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^ ; ABSTRACT 

. We investigate the physics driving the cosmic star formation (SF) history using the 

more than fifty large, cosmofogical, hydrodynamical simulations that together com- 
Q \ prise the Overwhelmingly Large Simulations (OWLS) project. We systematically vary 

■ the parameters of the model to determine which physical processes are dominant and 
' which aspects of the model are robust. Generically, we find that SF is limited by the 
, build-up of dark matter haloes at high redshift, reaches a broad maximum at inter- 
mediate redshift, then decreases as it is quenched by lower cooling rates in hotter and 
lower density gas, gas exhaustion, and self-regulated feedback from stars and black 
holes. The higher redshift SF is therefore mostly determined by the cosmological pa- 

• rameters and to a lesser extent by photo-heating from reionization. The location and 

Q>^ \ height of the peak in the SF history, and the steepness of the decline towards the 

,—1 . present, depend on the physics and implementation of stellar and black hole feedback. 

' Mass loss from intermediate-mass stars and metal-line cooling both boost the SF rate 

, at late times. Galaxies form stars in a self-regulated fashion at a rate controlled by 

■ the balance between, on the one hand, feedback from massive stars and black holes 
Q>^ ' and, on the other hand, gas cooling and accretion. Paradoxically, the SF rate is highly 

, insensitive to the assumed SF law. This can be understood in terms of self-regulation: 

i;^ ' if the SF efficiency is changed, then galaxies adjust their gas fractions so as to achieve 

the same rate of production of massive stars. Self-regulated feedback from accreting 
black holes is required to match the steep decline in the observed SF rate below red- 
\^ ' shift two, although more extreme feedback from SF, for example in the form of a 

5^ \ top-heavy initial stellar mass function at high gas pressures, can help. 

Key words: cosmology: theory - galaxies: formation - galaxies: evolution - stars: 
formation 
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1 INTRODUCTION 

The cosmic history of star formation (SF) is one of 
the most fundamental observables of our Universe. Mea- 
suring the global star formation rate (SFR) density 
as a function of redshift has therefore long been one 
of the pr i mary goals of observatio n al astronomy (e.g . 
iLillv et al.1 1 19961 : iMadau et all Il996l : ISteidel et all 1 19991 : 
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Ouchi et al. 2004 : Schim inovich et al. 20051: Arnouts et al.l 
20051 : iHopkins fc Beacom.200& : .Bouwens et al.ll200iD . 



Modeling the cosmic star formation history (SFH) is 
not an easy task because it depends on a complex inter- 
play of physical processes and because a large range of halo 
masses contributes. To predict the SFH within the context of 
the cold dark matter cosmology, one must first get the dark 
matter halo mass function right. These days this is the eas- 
ier part, as the cosmological parameters are relatively well 
constrained. Then one must model the rate at which gas 
accretes, cools, collapses, and turns into stars. Even if one 
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does not attempt to model the cold, interstellar gas phase 
and uses empirical SF laws to estimate the rate at which 
gas is converted into stars on kpc scales, there are a host of 
feedback processes that need to be taken into account. Stars 
produce radiation which can heat gas, exert radiation pres- 
sure, and change its ionization balance and hence its cooling 
rate. Massive stars explode as supernovae (SNe) which can 
drive both small-scale turbulence and large-scale outflows. 
Stars also produce heavy elements and dust which change 
the rate at which gas cools. Accretion onto supermassive 
black holes (BHs) in the centers of galaxies also results in 
radiative, thermal and mechanical feedback. Finally, mag- 
netic fields and cosmic rays may be important. 

Despite this complexity, many authors have 



used semi-analytic models (e.g. White & Frenk "l99ll: 



Hernguist fc Springel 20031: Bower et al..,2006.: iCroton et al.l 



20061 : ISomerville et al.1 l2008l : iFontanot et al.1 l2009l) 
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2002 


: Snrineel & HernauistI 2003bl: ISommer-Larsen et al. 
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NaEamine et al.l 


2009: 


OoDenheimer & Davel l2008l: 


Di Matteo et al.l 
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Choi & NaEamind 20091: ICrain et al.l 20091: Booth & Schave 


2009 


a|) to study the SFH of the Universe. 



Modern semi-analytic models use cosmological dark 
matter simulations to generate halo mass functions and 
merger trees. These are then combined with simple prescrip- 
tions for the baryonic physics, which include a large number 
of free parameters, to predict galaxy SFRs. The parameters 
of the models are tuned to match particular sets of obser- 
vations, after which the model is used to make predictions 
for other properties. While this approach has been very pro- 
ductive, it can be difficult to understand what physics is 
driving the results and there is a danger that model predic- 
tions may be correct even when the underlying prescriptions 
do not reflect the real world (because they were tuned to 
match particular observations). Furthermore, current semi- 
analytic models are only of limited use for the study of the 
gas around and in between galaxies. 

While hydrodynamical simulations attempt to model 
much more physics from first principles than semi-analytic 
models, they do make extensive use of subgrid prescrip- 
tions for the physics that they cannot resolve directly. Some 
of these prescriptions are physically motivated and rela- 
tively well understood. One example is the radiative cool- 
ing rate which is determined by atomic physics. However, 
even here there are significant uncertainties such as the ef- 
fects of photo-ionization, non-equilibrium, dust, and relative 
abundance variations. Other subgrid prescriptions, such as 
the stellar initial mass function (IMF) , are empirically moti- 
vated. Most subgrid models contain, however, a mixture of 
physical and empirical ingredients. For example, prescrip- 
tions for SF, thermal/kinetic feedback from stars and active 
galactic nuclei (AGN), and for the growth of BHs typically 
use physically motivated frameworks whose parameters are 
calibrated using observations. 

When it comes to predictions for SFRs in galaxies, the 
subgrid models used in cosmological hydro simulations play 
a critical role. Differences in the subgrid prescriptions used 
by different authors are likely to be much more important 
than differences in the codes used to model gravity and 
hydrodynamics. The crucial role played by subgrid mod- 



els implies that hydrodynamical simulations have some of 
the same weaknesses as semi-analytic models. On the other 
hand, direct simulation enables one to probe much deeper, 
e.g. by following gas flows in three dimensions, or by study- 
ing the interactions of galaxies with their environment and 
the intergalactic medium (IGM). 

While simulations offer the tantalizing possibility of 
laboratory-like control, the huge dynamic range required 
makes it computationally very expensive to explore parame- 
ter space. The Overwhelmingly Large Simulations (OWLS) 
project, which was made possible by the temporary avail- 
ability of the supercomputer that serves as the correlator 
for the LOFAR telescope, aims to use the potential of sim- 
ulations to gain insight into the physics that determines 
the formation of galaxies and the evolution of the IGM. 
The OWLS project consists of a large suite of cosmolog- 
ical, smoothed particle hydrodynamics (SPH) simulations 
with varying box sizes and resolutions. Each production sim- 
ulation uses 2 X 512'' particles, which places them among 
the largest dissipative simulation ever performed. The real 
power of the project stems, however, not from the size of the 
simulations, but from the fact that they are repeated many 
(more than 50) times, each time varying a subgrid prescrip- 
tion, most of which were newly developed for this project. 
Although we have not attempted to fine-tune the subgrid 
parameters to match particular data sets, we do hope that 
our investigations will help trigger future work in directions 
that will improve agreement with the observations. 

While the OWLS project aims to study a variety of 
problems in astrophysical cosmology, we will limit ourselves 
here to the cosmic SFH. To further limit the scope of the 
paper, we will postpone a discussion of the SFR as a function 
of halo mass to another paper (Haas et al., in preparation). 
We note, however, that this latter function is in some ways 
physically more interesting than the SFH. This is because 
the cosmic SFR can be thought of as a convolution of the 
halo mass function, which depends only on cosmology and 
redshift, and the SFR as a function of halo mass and redshift, 
which depends on poorly understood astrophysics. 

We will also not compare to observations of the build 
up of the cosmic stellar mass density. The cumulative SFR 
must of course equal the stellar mass (after taking stellar 
mass loss into account). While this will be true by construc- 
tion for ab initio models, it is not necessarily the case for 
observationally inferred quantities. This is because observa- 
tional probes of SF measure only the rate of formation of 
massive stars, since those dominate the electromagnetic out- 
put of young stellar populations. To obtain the total SFR, it 
is necessary to extrapolate to low masses using an assumed 
IMF. Since the IMF is uncertain and may not be universal, 
measurements of the total stellar mass are of great interest. 
However, it should be noted that massive stars are the ones 
that matter from a cosmological perspective because they 
dominate the chemical, radiative and mechanical feedback 
processes that regulate the formation of stars and galaxies. 
As the objective of this paper is to explore the effects of 
varying physical prescriptions rather than to fit the obser- 
vations, we will for simplicity limit ourselves to comparisons 
with measurements of the evolution of the SFR density. We 
emphasize, however, that it should be kept in mind that 
these are often in consistent with observations of the stellar 
mass density (e.g. lWilkins et al ]|2008l : ICowie fc Bar"^l2008l . 
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Table 1. List of simulations run assuming the reference model. From left-to-right the columns show: simulation identifier; comoving box 
size; number of dark matter particles (there are equally many baryonic particles); baryonic particle mass; dark matter particle mass; 
comoving (Plummer-equivalent) gravitational softening; maximum physical softening; final redshift. 



Simulation 


L 

(/i-i Mpc) 


AT 


(h-iMo) 


m-dm 


Ccom 

(/i-i kpc) 


£prop 

(h-^kpc) 


^cnd 


REF.L006N128 


6.25 


128^ 


1.4 X 10*^ 


6.3 X lO*' 


1.95 


0.50 


2 


REF.L012N256 


12.50 


256^ 


1.4 X 10*^ 


6.3 X 10^ 


1.95 


0.50 


2 


REF.L012N512 


12.50 


512^ 


1.7 X 10^ 


7.9 X 10^ 


0.98 


0.25 


2 


REF.L025N128 


25.00 


128^ 


8.7 X 10^ 


4.1 X 10* 


7.81 


2.00 





REF.L025N256 


25.00 


256=^ 


1.1 X 10^ 


5.1 X lO'^ 


3.91 


1.00 


2 


REF.L025N512 


25.00 


512^ 


1.4 X lO'' 


6.3 X 10^ 


1.95 


0.50 


1.45 


REF.L050N256 


50.00 


256=^ 


8.7 X 10^ 


4.1 X 10* 


7.81 


2.00 





REF.L050N512 


50.00 


512^ 


1.1 X 10^ 


5.1 X lO'^ 


3.91 


1.00 





REF.L100N128 


100.00 


128^ 


5.5 X 10^ 


2.6 X 10^0 


31.25 


8.00 





REF.L100N256 


100.00 


256^ 


6.9 X 10* 


3.2 X 10^ 


15.62 


4.00 





REF.L100N512 


100.00 


512^ 


8.7 X 10^ 


4.1 X 10* 


7.81 


2.00 






but see also iReddv at al1l2008l : iReddv fc Steidelll2009l) and 
subje ct to large systematic uncertainties (e.g. IConrov et al.l 
I2OO9I ). 

This paper is organized as follows. In section [2] we in- 
troduce the numerical techniques that are common to all 
OWLS runs, while sections|3]and|4]describe the physics mod- 
ules employed in the reference model and the other OWLS 
runs, respectively. As this paper also serves to introduce the 
OWLS project, we will discuss the ingredients of the differ- 
ent simulations in some detail. The SF histories predicted 
by the models are presented and discussed in the sections 
that describe them. Finally, we discuss and summarise our 
main findings in section [S] 



2 SIMULATIONS 

In this section we describe the numerical techniques that are 
common to all simulations that make up the OWLS project. 
The modifications and additions to the subgrid modules are 
discussed in section|3]for our reference model and in section|4] 
for the other runs. 

The simulations were performed with a significantly ex- 
tended version of the i V-Body Tree-P M, SPH code gad- 
gets (last described in ISpringell |2005| ) , a Lagrangian code 
used to calculate gravitational and hydrodynamical forces 
on a system of particles. Most simulations were run in pe- 
riodic boxes of size L = 25 and 100 comoving Mpc 
and each of the production runs uses 512^ dark matter and 
equally many baryonic particles (representing either colli- 
sionless star or collisional gas particles). The particle masses 
in the 2 x 512^ particle 25 Mpc (100 h^'^ Mpc) box are 

6.34 x lO** h-^ Mq (4.06 x 10* h'^ Mq) for dark matter and 

1.35 x W^h-^MQ (8.66 x IO'^^-^Mq) for baryons. Note, 
however, that baryonic particle masses change during the 
course of the simulation due to mass transfer from star to 
gas particles. The 25 Mpc simulation volumes were run 
as far as redshift 2 and the 100 Mpc volumes were run 
to redshift 0. Comoving gravitational softenings were set to 
1/25 of the initial mean inter-particle spacing but were lim- 
ited to a maximum physical scale of 0.5 kpc/h (2 kpc/h) 
for the 25 Mpc (100 Mpc) simulations. The switch 
from a fixed comoving to a fixed proper softening happens 



at z = 2.91 in all simulations. We used A'ngb = 48 neighbors 
for the SPH interpolation. 

In order to assess the effects of the finite resolution and 
box size on our results, we have run a suite of cosmological 
simulations, all using the same physical model, using differ- 
ent box sizes (ranging from 6.25 Mpc to 100 Mpc) 
and particle numbers (ranging from 128^ to 512^). The par- 
ticle masses and gravitational softenings for each of these 
simulations are listed in Table [T] 

The ini tial conditions were gener ated with CMBFAST 
(version 4.1; ISeliak fc Zaldarriagal 19961 ) and evolved to red- 
shift z = 127, whe re the simulations were started, using the 
IZel'Dqvichl (Il970l) approximation from an initial glass-like 
state (|Whitdll996l ). 

Before discussing each of the variations made to the 
subgrid models in section |3J we first turn out attention to 
the physics included in the reference model. 



3 THE REFERENCE MODEL 
3.1 Description of the model 

When simulations lack the required numerical resolution or 
the physics to accurately model a physical process, we must 
resort to subgrid models. In this section we describe the 
'reference' physical model (REF) that is used as a base 
from which all further investigations of the behaviour of the 
simulations can be launched. As discussed in section [J] we 
will do this by varying one physical process at a time and 
comparing the resulting SFH to the one predicted by the 
reference model. 

We emphasize that the REF model should not be 
viewed as our "best" model. As its name implies, it func- 
tions as a reference point for our systematic variation of 
the parameters. In fact, we will show in future papers that 
some variations yield much better agreement with particu- 
lar types of observations. For example, the inclusion of AGN 
feedback dramatically improves the agreement with obser- 
vations of groups of galaxies at redshift zero (McCarthy et 
al., in preparation). 

To illustrate the dynamic range in the simulations. 
Figs. [T] and [5] show two factors of ten zooms into the tenth 
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most massive haloes in, respectively, the REF_L025N512 
run at z = 2 and model REF.L100N512 aX z = Q. Note, 
however, that the left panels still show only a small fraction 
of the simulation volume. 



3.1.1 Cosmology 

We assume the values for the cosmological parame- 
ters derived from the Wilkinson Microwave An is otrop y 
Probe (WMAP) 3-year results (|Spergel et al.1 |2003), 
{nm,nb,r2A,a8,na,/i} = {0.238, 0.0418, 0.762, 0.74, 0.951, 
0.73}, which are con sistent with the WMAP 5-year data 
jKomatsu et al.l|200 8f) Fl The primordial baryonic mass frac- 
tion of helium is assumed to be 0.248. 



3.1.2 Radiative cooling and heating 

Radiative cooling is central to simulations of the formation 
of galaxies as it enables baryons to dissipate their binding 
energy which allows their collapse to proceed within virial- 
ized structures. Photo-heating by the ionizing background 
radiation also plays a key role because it strongly increases 
pressure forces in low-density gas, thereby smoothing out 
small-scale baryonic structures. 

Previous cosmological simulations have typically in- 
clude d radiative cooling assuming primordial abundances 
(e.g. ISpringel fc Hernguistl l2003bl '). So me recent studies 
have ir icluded meta l- line c o oling (e.g. Scannapieco et all 



20051: iRomeo et al.1 l2006l: [Op ponhcimcr fc Pavel 1200^ 



Tornatore et al.ll2007l : IChoi fc Na ga niine ,2009 '). but under 



the assumption of coUisional ionization equilibrium and 
fixed relative abundances. Photo-ionization by the UV back- 
ground radiation does not only provide a source of heat, it 
also reduces the cooling rates for bo t h primordial and meta l- 
enriched plasm a s (Efst athioul 1 19921 : IWiersma et al.l l2009al ) . 
|w iersma et al.l (|2009al ) emphasized the importance of in- 
cluding this effect as well as variations in the relative abun- 
dances of the elements. 

We implemented radia tive cooling and heating using the 
method and tables of Wiersma et aLri|2009a 'FI. In brief, net 
radiative cooling rates are computed element-by-element in 
the p resence of the cosmic microwave background and the 
iHaardt fc Madau (2001) model for the UV and X-ray back- 
ground radiation from quasars and galaxies. Hence, varia- 
tions in relative abundances and photo-ionization of heavy 
elements are both taken into account. The contributions of 
the eleven elements hydrogen, helium, carbon, nitrogen, oxy- 
gen, neon, magnesium, silicon, sulphur, calcium, and iron 
are interpolated as a function of density, temperature, and 
redshift from tables that have been precomputed using the 
publicly avai lable photo-ionization package CLOUDY, last 
described bv lFerland et al.] l| 19981 ). assuming the gas to be 
optically thin and in (photo-)ionization equilibrium. 

The simulati ons model hydrogen reionization by 
'switching on' the iHaardt fc Madaul (|200ll ) background at 
2 = 9. Prior to reionization the cooling rates are computed 



^ The most notable difference is in erg, which is 1.6 cr lower in 
WMAP3 than in WMAP5. 

^ We used their equation (3) rather than (4) and Cloudy version 
05.07 rather than 07.02. 



in the presence of the cosmic microwave background and 
a photo-dissociating background whic h we o btain by cut- 
ting off the z = 9 IHaardt fc Madaul (|200ll ) spectrum at 
1 Ryd. Note that the presence of a photo-dissociating back- 
ground suppresses H2 cooling at all redshifts. Reionization 
has the effect of rapidly heating all of the gas to tempera- 
tures ~ 10"* K. The assumption that the gas is optically thin 
is likely to lead to an underest imate of the gas temper ature 
shortly after reionization (e.g. lAbel fc Haehneltlll999l ). For 
the case of helium reionization we correct for this effect by 
heating the gas by a total amount of 2 eV per atom. This 
extra helium reionization heating takes place at a central 
redshift of 3.5, with the heating spread with a Gaussian filter 
with (j{z) = 0.5 in redshift. This prescription was chosen to 
match observations of the temperature hist ory of the IGM 
' Schave et all I2OO0I ) as shown in Fig. 1 of IWiersma et al.l 



2009bl ). 



3.1.3 Star formation 

Cosmological simulations such as ours miss both the resolu- 
tion and the physics to model the cold interstellar medium 
(ISM), let alone the formation of stars within molecular 
clouds. Star formation is therefore implemented stochasti- 
cally by converting gas particles into coUisionless star parti- 
cles, which represent simple (or single) stellar populations. 
We convert entire particles, because the spawning of multi- 
ple star particles p er gas particle affects the effi ciency of 
feedback from SF (|Dalla Vecchia fc Schavl I2OO8I ). Hence, 
the particle number is conserved in our simulations. 

Gas with densities exceeding the critical density for the 
onset of the thermo- gravitational instability (hydrogen num- 
ber densities riH = 10~^ — lO "'^ cm~ ' ^) is expected to be 
multiphase and to form stars l|Schavell2004h . We therefore 
impose an effective equation of state (EOS) with pressure 
P (X p^"*' for densitieijfl nn > J^H where rin = 0.1 cm~'^, 
normalised to P/k = 1.08 x 10^ cm"^ K at the threshold. We 
use 7cff = 4/3, for which both the Jeans mass and the ra- 
tio of the Jeans length to the SPH kernel are independent of 
the density, thus preventing s purious fragmentation due to a 
lack of numerical resolution jSchave fc Dalla Vecchia|[2008l : 
iRobertson fc Kravtsovll2008l). Only gas on t he EO S is al- 
lowed to form stars. Schave fc Dalla Vecchial (|2008l ) demon- 
strated that our choice of threshold reproduces the threshold 
surface density for Ha emission from SF that is observed in 
nearby galaxies. 

Previous cosmological simulations used Schmidt-type 
(i.e., power-laws of the volume density) SF laws and tuned 
one or more free parameters to fit the observed Kennicutt- 
Schmidt SF law, which is a surface density law. This ap- 
proach is unsatisfactory as the parameters would really need 
to be re-tuned if disk scale heights change as a result of vary- 
ing abundances (and hence cooling rates), SN feedback or 
cha nges in the assumed EOS of th e ISM. 

ISchave fc Dalla Vecchial l|2008l ) showed that because the 



^ Gas particles are only placed on the EOS if their temperature 
was below 10^ K when they crossed the density threshold and if 
their density exceeds 57.7 times the cosmic mean. These criteria 
prevent SF in intraclu ster gas and in intergalactic gas at very high 
redshift, respectively l lSchave &: Dalla Vecchiall200a) . 
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Figure 1. Zoom into a M200 = 10^^'^ Mq halo at z = 2 in tiie REF.L025N512 simulation. Prom left-to-right, the images are 10, 1, and 
0.1 Mpc on a side. All slices are 1 Mpc thick. Note that the first image shows only a fraction of the total simulation volume, 
which is cubic and 25 Mpc on a side. The color coding shows the projected gas density, logjQ p/ (pb)i a^nd the color scale ranges from 
-1 to 4 (which is lower than the true maximum of the image). The coordinate axes were rotated to show the galaxy face-on. This halo 
is the 10th most massive in the simulation. About half of the haloes in this mass range host extended disk galaxies, while the other half 
have highly disturbed morphologies due to ongoing mergers. 




Figure 2. Zoom into a M200 = lO^'''^ Mq halo at 2 = in the REF.L100N512 simulation. From left-to-right, the images are 40, 4, and 
0.4 Mpc on a side. All slices are 1 Mpc thick. Note that the first image shows only a fraction of the total simulation volume, 
which is cubic and 100 Mpc on a side. The color coding shows the projected gas density, logj^Q p/ (pb)i and the color scale ranges 
from -1 to 4 (which is lower than the true maximum of the image). This halo is the 10th most massive in the simulation. 



surface density in a self-gravitating system is directly re- 
lated to the pressure, the Kennicutt-Schmidt SF law can 
be rewritten as a pressure law. This enables one to re- 
produce arbitrary input Kennicutt-Schmidt laws indepen- 
dently of the assumed EOS. Moreover, because the pa- 
rameters are observable s , no t uning is required. Following 
ISchave fc Dalla Vecchial l|2008l ). we thus compute the SFR 
for star-forming gas particles using 



,^(lMe pc-^)-"(g/,p) 



(n-l)/2 



(1) 



where rrig is the mass of the gas particle, 7 = 5/3 is the 
ratio of specific heats (not to be confused with the effec- 
tive EOS imposed onto the ISM), fg is the mass fraction in 
gas (which we assume to be unity) and P is the total pres- 
sure. The parameters A and n are, respectively, the ampli- 
tude and slope of the observed iKennicuttI ()l998l ) law, E* = 



A(Eg/l Mq pc"^)" with A = 1.515 x 10"'' Mgyr-^kpc"^ 
and n = 1.4. The amplitude of this relation has been renor- 
malised by a factoi FI 1/1.65 to accou nt for the fact that the 
origin al analysis of KennicuttI lll998l) assumed t he ISalpeteJ 
(jl955h IMF whereas we use the lChabrie^ l|2003l ) IMF. 



3.1.4 Stellar evolution and chemodynamics 

Our implementation of stellar evo lution and chemical en - 
richment is discussed in detail in lersma et al.l (|2009bh . 
Here, we will provide only a brief summary. 

Each star particle represents a single stellar population 



* This normalization factor is calculated from the asymptotic ra- 
tio (which is reached after only 10* yr) of the numbers of ionizing 
photons pr edicted from models of ste llar populations with a con- 
stant SFR l lBruzual fc Charlotibooal') . 
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that is specified by its initial mass, age, and its chemical 
composition (which it inherits from its progenitor gas par- 
ticle). We follow the timed release, by both massive stars 
(Type II SNe and stellar winds) and intermediate mass stars 
(Type la SNe and asymptotic giant branch (AGB) stars), of 
all 11 elements that contribute significantly to the radiative 
cooling rates. During each time-step, star particles distribute 
the mass they eject over their neighboring gas particle^ us- 
ing the SPH inter polation scheme . 

We assume a Ichabriej (|2003| ) IMF spanning the range 
0.1 to 100 Mq and use t h e me tallicity-dependent stellar 
lifetimes of iPortinari et aL (Il998l) and the complete set of 
nucleo synthetic yields of Marigol 1 200lh and IPortinari et al.l 
l| 19981 ) al ong with the SN Type la (SNIa) yields of the W7 
model of iThielemann et al.1 (120031 ). Since SNIa are thought 
to result from binary evolution, a single stellar population 
will produce SNIa over an extended period. We implement 
the release of mass and energy (which we inject in thermal 
form) by SNIa using empirically derived rate s normalised to 
the ob served cosmic SNIa rate (see Fig. A6 o f|w iersma et al.l 
l2009bl ). For reference, our assumed rate implies that a frac- 
tion of about 0.025 of stars with initial mass between 3 and 
8 solar masses end their lives as Type la SNe. 

For the purpose of both radiative cooling and stellar 
evolution, we define the abundance of a particular element 
as the ratio of the SPH estimates of its mass de nsity and the 
total gas mass density, iw lersma et al.l (|2009bh showed that 
the use of such 'smoothed abundances' for the cooling rates 
significantly increases the SFR compared to the standard 
approach of using 'particle metallicities', i.e., the ratio of 
the elemental mass to the total mass of a particle. While 
the use of smoothed abundances reduces the effects of the 
lack of metal mixing inherent to SPH, it does not solve the 
problem. 



3.1.5 Energy feedback from core collapse supernovae 

The use of an effective EOS with a pressure that exceeds 
that of the warm, neutral ISM can be considered a form of 
weak feedback that reflects the fact that energy injected by 
massive stars and SNe drives small-scale turbulence. How- 
ever, as we will show explicitly in section 14.41 this form of 
feedback does not lead to a signiflcant suppression of SF. In- 
deed, observatio ns show th at starburst galaxies drive large- 
scale winds (e.g. iVeilleux e t al. 2005) which may, over time, 
eject large amounts of gas and may therefore dramatically 

reduce the SFR. 

As discussed in lDalla Vecchia fc Schav3 (|2008l ). thermal 
energy from SNe is quickly radiated away in simulations like 
ours because the ratio of t he heated mass to that of the 
star particle is too large. In iDalla Vecchia fc Schavj (|2009l ) 
and section 14.8.31 we show that it is possible to overcome 
this overcooling problem, without ad-hoc suppression of the 
radiative cooling rates, by decreasing this ratio. However, for 
our reference model we use the more standard approach of 
injecting SN energy in k i netic f orm using the prescription of 
iDalla Vecchia fc Schav3 (1200 8'). which is a variation of the 
recipe of ISpringel fc Hernguist (2003a ). 



^ As discussed in lWiersma et aL I ll2009d) . we do not change the 
entropy of the receiving particles. 



After a short delay of 30 Myr, corresponding to the max- 
imum lifetime of stars that end their lives as core-collapse 
SNe, newly-formed star particles inject kinetic energy into 
their surroundings by kicking a fraction of their SPH neigh- 
bours in random directions. Wind particles are not allowed 
to form stars for a period of 15 Myr in order to avoid high 
velocity star particle ejection (a numerical artifact we ob- 
served occasionally in high-resolution simulations of isolated 
galaxies). These time delays are not important for the results 
presented here. 

Each SPH neighbour i of a newly-formed star particle j 
has a probability of rynijl '^}2n-i' receiving a kick with 

a velocity Thus, if all baryonic particles had equal mass, 
each newly formed star particle would kick, on average, t] 
of its neigh bours. Our reference model use s the default pa- 
rameters of iDalla Vecchia fc Schav3 (|2008l ). i.e. r\ = 1 and 
Ww ~ 600 kms"'^ Assuming that each star with initial mass 
in the range 6 — 100 M© injects 10^^ erg of kinetic energy, 
these parameter values imply that the total wind energy ac- 
counts for 40 per cent of the available kinetic energy for our 
IMF (if we ignore the electron capture SNe predicted by 
models with convective overshoot (e.g. IChiosi et al.l Il992l ) 
and consider only stars in the mass range 8 — 100 M© , this 
works out to be 60 per cent). The value r\ = 2 was chosen in 
part because it roughly reproduces the peak in the cosmic 
SFR, as we will show later. 

Note that contrary to the widely- u sed kinetic feedback 
recipe of ISpringel fc Hernguistl (l2003al ). the kinetic energy 
is injected locally and the wind pa rticles are not decoupled 
hydro dynamically. As discussed bv lDalla Vecchia fc Schav3 
(|2008l ) and as we will show in Section |4]8]2] these differences 
have important consequences. 



3.1.6 Black hole growth and feedback from AGN 

The reference model does not include a prescription for the 
growth of supermassive BHs and feedback from AGN. How- 
ever, AGN are included in the OWLS project as an optional 
model that is switched on for a subset of the simulations 

f Sec, mm. 



3.2 Convergence tests 

Before trying to interpret the results of numerical simu- 
lations, we must check whether they have converged nu- 
merically. For cosmological simulations this means checking 
whether the box size was sufficiently large and whether the 
resolution was sufficiently high. To isolate the effects of the 
size of the simulation volume and the resolution, it is neces- 
sary to vary one while holding the other fixed. 

To test for numerical convergence we have run a suite 
of simulations of the reference model. The main numerical 
parameters of these runs are listed in Table [T] The simula- 
tion names contain strings of the form LxxxNyyy, where xxx 
is the simulation box size in comoving Mpc and yyy is 
the cube root of the number of particles per species (dark 
matter or baryonic). 

For readers not interested in the details, we first give 
the conclusions so that they can skip to section |4] Even 
our 25 Mpc boxes are sufficiently large to obtain a con- 
verged prediction for the cosmic SFH. The resolution of the 
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Figure 3. The effect of the box size on the cosmic SFH. The 
curves show the cosmic SFR density as a function of rodshift 
(bottom X-axis) and lookback time (top x-axis) for different sim- 
ulations of the reference m odel. The data point s show the compi- 
lation of observations from lHopkins fc Beacoml | |2006| ). converted 
to our IMF and cosmology. The black curves extending down 
to z = are for simulations that use the same numerical reso- 
lution as model REF.L100N512. The red curves, which do not 
continue to z = 0, correspond to runs with the resolution of 
REF.L025N512. A 25 h'^ Mpc box is sufficiently large to ob- 
tain a converged prediction for the SFH down to 2 = 0. 

L025N512 and L100N512 runs suffices for redshifts z <7 
and z < 3, respectively. While the SFR typically increases 
if the resolution is improved, the situation reverses at low 
redshift before convergence has been attained. 

The disc ussion below is i n pa rts similar to the one we 
presented in IWiersma et al. 

I (|2009b), where we considered 
the convergence of the cosmic metal distribution in the ref- 
erence simulations. 

3.2.1 Box size 

The results of cosmological simulations may depend on the 
size of the simulation volume for at least two reasons. First, 
because the mean density in the box is fixed by the cosmo- 
logical parameters, the box size determines what types of 
objects can be sampled. The distribution of density fluctua- 
tions can only be modeled correctly on scales that are much 
smaller than the box. Second, because the Fourier modes 
of the density field only evolve independently in the linear 
regime, the box must be large compared to the scales on 
which the density contrast is non-linear. Otherwise the miss- 
ing power on scales greater than the box size will decrease 
the power on the scales that are sampled by the simulation. 

Fig. |3] shows the SFR per unit comoving volume 
as a function of redshift for two sets of simulations 
of the reference model. The solid curves show our two 
fiducial box sizes and particle numbers: 100 comoving 
/i"^ Mpc (REF.L100N512, black) and 25 comoving Mpc 
{REF.L025N512, red), both using 2 x 512^ particles. 
The data points s h ow t he observations as compiled by 
iHopkins fc Beacoml (|2006D . To facilitate easy comparisons, 
we will show these data points and at least one of the two 
fiducial runs in all subsequent figures. 

We caution the reader that the data are subject to large 



systematic uncertainties due to for example the assumed 
IMF and dust correction. Observe also that the scatter is 
clearly t oo large compared w i th th e error bars, despite the 
fact that iHopkins fc Beacoml (|2006l ) applied a uniform dust 
correction and that the same IMF was assumed for all obser- 
vations. Given these uncertainties, models whose predictions 
are discrepant with respect to these observations cannot au- 
tomatically be ruled out. 

Focusing first on the LlOO run (black, solid curve) we 
see a sharp rise at high redshift, a peak at z ~ 2 followed 
by a steady decline to z = 0. Qualitatively this matches 
the data, although the simulations appear to underestimate 
the SFR beyond the peak as well as the steepness of the 
decline below 2 = 1. Note that the height of the peak, i.e. 
the maximum SFR, is sensitive to the fraction of the SN en- 
ergy that is used to generate galactic winds. While we fixed 
the wind velocity to 600 kms~^ based on other considera- 
tions (see iDalla Vecchia fc Schavell200^ . the mass loading 
rj — 2 (which corresponds to 40 percent of the SN energy 
for Uw = 600 kms~^), was chosen partly because it gives 
roughly the right maximum SFR. We will show later that 
the underestimate of the SFR at z > 3 can be attributed 
to a lack of numerical resolution, while the overestimate of 
the SFR at z < 0.5 refiects the fact that our galactic winds 
cannot suppress SF in massive galaxies. 

Comparing the three black curves extending to z = 0, 
which correspond to box sizes of 25, 50, and 100 h'^ Mpc, we 
see that even a 25 Mpc box is large enough to obtain a 
converged estimate for the cosmic SFH. This is perhaps sur- 
prising, as there clearly exist structures with sizes that are 
of the same order or greater than this. Apparently, rare ob- 
jects like clusters of galaxies do not con tribute significantl y 
to the mean SFR. This is consistent with lCrain et al.| l|2009l ). 
who used zoomed simulations to show that while the SFR in 
different 25 Mpc regions varies by up to an order of mag- 
nitude, the SFH in a region of this size whose mean density 
equals the cosmic mean closely tracks the global SFH. 

Comparing the three red curves that end at higher red- 
shifts, which correspond to box sizes of 6.125, 12.5 and 
25 Mpc and particle masses that are 64 times smaller 
than those used for the black curves, we see that while a 
12.5 Mpc box is nearly sufficiently large for z > 2 (recall 
that we already established that 25 Mpc is sufficiently 
large using the lower resolution simulations), 6.125 Mpc 
is clearly insufficient to obtain a converged estimate of the 
SFH. 

Now that we have established that our fiducial box sizes 
are sufficiently large, we turn our attention to the conver- 
gence with respect to resolution. 

3.2.2 Numerical resolution 

While the simulation box limits the maximum sizes of the 
structures that can form, numerical resolution may even af- 
fect the properties of common objects. For example, hydro- 
dynamical simulations that do not resolve the Jeans scales 
may underestimate the fraction of mass in collapsed struc- 
tures and thus the SFR. Moreover, we cannot expect to 
form dark matter haloes whose masses are comparable to 
or smaller than the particle mass. 

Before showing the results of the convergence tests, it 
is useful to consider what to expect. The L025N512 model 
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Figure 4. As Fig. |3] but comparing the SFHs for simulations of the reference model that differ in terms of their numerical resolution. 
The left and right panels test the convergence of the SFH predicted by the L025N512 and L100N512 runs, respectively. The SFH in 
the L025N512 run is clearly fully converged for z < 4 (compare with the lower resolution L025N256), but comparison with the higher 
resolution L012N512 (which, however, has a box size that is slightly too small to be converged; see Fig. [Sj suggests that it is nearly 
converged for z < 7. Comparing L100N512 with L050N512 (which uses a box size sufficiently large to provide a converged result; see 
Fig-EJi W6 see that the former has nearly converged for z < 3. Interestingly, increasing the resolution beyond that of L100N256 decreases 
the SFR at late times. 



has a dark matter particle mass of mdm ~ 6 x 10^ Mq. 
From a comparison of the mass functions of dark matter 
only simulations, we find that ~ 10^ particles are needed 
to robustly define a halo. Thus, we expect the halo mass 
function to be converged for M > 6 x 10^h~^MQ. Com- 
paring this to the mass corresponding to the virial tempera- 
ture below which photo-heating is expected to suppress SF, 
M ~ 0.7 X 10« h-^ Mq (Tvi./lO'* K)^''^ ((1 + ^)/10)-='/^ we 
see that even after reionization we are missing haloes ex- 
pected to form stars. However, haloes of such low mass are 
only expected to dominate the cosmic SFR at very high red- 
shift. The particle mass for L100N512 is 64 times greater 
and we can thus only probe the mass function down to 
M ~ 4 X 10^" /i"^ Mq. This wiU l ead us to un derestimate 
the SFR at 2>3 (see e.g. Fig. 5 of lCrain et aL. 2009). 

Let us now consider the hydrodynamics. The Jeans 
scales depend on the density and the temperature of the gas. 
The temperature of substantially overdens^fl gas is > 10* K 
in our simulations. In reality, gas at interstellar densities 
(riH ^ 10~^ cm~^) is sufficiently dense and self-shielded to 
form a cold (T < 10'' K) gas phase, allowing it to form stars 
l|Schavdl2004l ). However, our simulations impose an effective 
EOS for gas with densities exceeding our SF threshold of 
riH = 10~'^cm~^. For our EOS (P oc p*^^) the Jeans mass 
is independent of the density. Hence, if we resolve the Jeans 
mass at the SF threshold, we resolve it everywhere. The 
Jeans mass is given by 

(tF?I^) (k?k) 

where /g is the (local) fraction of the mass in gas. Thus, 

^ Gas with very low overdensities can have temperatures substan- 
tially below 10* K due to the adiabatic Hubble expansion, but the 
Jeans scales corresponding to these low densities are nevertheless 
large. 



we do not expect convergence unless the gas particle mass 
rUg <^ 10^ Mq. To achieve convergence, a simulation will, 
however, also need to resolve the Jeans length Lj. This 
implies that the maximum, proper gravitational softening, 
Eprop, must be small compared with the Jeans Length 

'^-/.""(ra)"""(llK)"'' 

Note that since Lj scales as Lj oc p~^^^ for our EOS, eprop 
will always exceed Lj for sufficiently high densities. How- 
ever, since the Jeans mass does not decrease with density, 
we do not expect runaway collapse for star-forming gas. 

Comparing the above equations with the gas particle 
mass and softening scales for our fiducial simulations (see 
Table [l|, we see that while L025N512 marginally resolves 
the Jeans scales for /g « 1, this is not the case for the sim- 
ulations that go down to 2 = 0, although L050N512 has 
mg « Mj and tprop ~ L} and is therefore not far off. Note, 
however, that none of our simulations come close to resolv- 
ing the Jeans scales prior to reionization, when SF in haloes 
with virial temperatures less than 10* K may have been im- 
portant. 

Fig. |4]compares the SFHs predicted by simulations with 
varying resolutions. The left and right panels test the con- 
vergence of the L025N512 and L100N512 models, respec- 
tively. Focusing first on the solid and dashed black curves 
in the left panel, we see that a particle mass 8 times greater 
than our fiducial value (and a softening twice our fiducial 
value) is sufficient for 2 < 4. For 2 < 2 even a particle 
mass that is a factor 64 smaller appears to be sufficient. 
Comparison of our fiducial run with the higher resolution 
run L012N512 (red, dot-dashed), indicates that the former 
is likely nearly converged for 2 < 7. Note that we do not 
expect perfect agreement when comparing L025 and L012 
runs even if they have converged in terms of resolution since 
the two runs necessarily have different initial conditions and 
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a 12.5 h Mpc box is not fully converged in terms of box 
size (see Fig. [3]). 

The kink at « = 9 in the SFH of L012N512 is due to the 
negative feedback associated with reheating during reioniza- 
tion. While this effect is easily visible for L012N512, it is 
only just detectable in L025N512. This is expected. Reion- 
ization will suppress SF in haloes with virial temperatures 
< 10* K, which corresponds to halo masses of ~ 10* M© at 
this redshift. Such haloes are resolved with ~ 10^ particles 
in L012N512, but contain only ~ 10 particles in L025N512. 

Comparing the three black curves in the right panel, 
which show the SFH for runs L WON 128, L WON 25 6, and 
LW0N512, we see no evidence for full convergence, although 
the difference between the two highest resolution runs is 
small for z < 1.5. Indeed, comparison with L050N512 and 
L025N512, whose box sizes we have already shown to be 
sufficiently large, reveals that the fiducial run LW0N512 has 
nearly converged for z < 3 (and L050N512 for 

An increase in the numerical resolution typically in- 
creases the SFR, particularly at high redshift, when it is 
dominated by haloes near the resolution limit. Observe, how- 
ever, that the opposite happens at low redshift once the reso- 
lution is increased beyond that of LW0N256. This reduction 
reflects the fact that, in simulations with higher resolution, 
the gas that would otherwise be available or SF has already 
been used up or ejected by lower mass progenitors at higher 
redshifts. 

Single simulations currently lack the dynamic range to 
obtain a converged result for the SFH over a wide range of 
redshifts. One s trategy that has been used to overcome this 
limitation (e.g. ISpringel fc Hernauist|[2003bl ) is to combine 
a suite of simulations with different box sizes and to use, 
for each redshift, the one that yields the highest SFR. This 
procedure comes down to plotting all SFH curves and using 
the envelope that encompasses all as the best estimate for 
the SFH. However, our finding that, at low redshift, the pre- 
dicted SFR decreases as it approaches convergence, indicates 
that this procedure may overestimate the SFR at late times. 
Because the low-redshift SFR depends on the amount of gas 
that was consumed or ejected at earlier times, one really 
does need a large dynamic range to model the SFH down to 
z = 0. 

Summarizing, the convergence tests are consistent with 
our expectations based on our estimates of the minimum 
resolved dark halo mass and a comparison of the mass and 
length resolutions with the Jeans scales. For our purposes, 
the resolution of the L WON 5 12 run suffices for 2: < 3 and 
that of L025N512 for z<7. 

Before investigating the physics driving the predicted 
SFH, we caution the reader that convergence of our refer- 
ence model does not automatically imply that other physics 
variations are also converged at the same resolution. On the 
other hand, it would be surprising if resolution effects would 
change qualitative conclusions drawn from comparisons in 
the regime for which REF is converged. This situation would 
change, however, if we did not impose an effective equation 
of state onto the ISM because in that case the Jeans scales 
could become much smaller than in the models considered 
here. 



4 VARIATIONS ON THE REFERENCE 
MODEL 

In this section we describe the full set of OWLS runs. Most 
simulations differ from the reference model in only the choice 
of a single parameter or the presence of a certain aspect of 
the subgrid physics. In this way, cross-comparison of dif- 
ferent simulations allows us to isolate the effects of different 
physical processes and the importance of different numerical 
parameters. The full list of simulations is shown in Table |2l 
which lists the simulation identifier, indicates whether or not 
a given simulation was run in the 25 and 100 Mpc boxes 
and gives a reference to the section that discusses the simu- 
lation. Except for the MILL runs, all simulations that were 
run in the same box size used identical initial conditions. 

The order in which we present the different variations 
roughly parallels the order in which the subgrid models were 
discussed in section 13.11 The different subsections can be 
read independently of each other. We begin by comparing 
our WMAP-3 cosmology to that of the WMAP-1 cosmology, 
which was for exa mple assumed in the w idely used "Millen- 
nium simulation" ijSpringel et al ]|2005b'). In sections l4.2l and 
I4.3l we investigate aspects of the radiative cooling and heat- 
ing by turning off metal-line cooling and varying the redshift 
of reionization, respectively. We study the importance of our 
treatment of the unresolved ISM in section 14.41 by varying 
the EOS. We vary the subgrid model for SF in section B31 
where we try a metallicity-dependent SF threshold as well 
as a range of Kennicutt-Schmidt SF laws. In section |4^ we 
investigate the effect of intermediate mass stars by turning 
off mass loss from AGB stars and by varying the time de- 
lay function for Type la SNe. Section 14.71 investigates the 
effect of using a Salpeter IMF and an IMF that becomes 
top-heavy at high pressures. Many aspects of our prescrip- 
tion for kinetic SN feedback are varied in section 14.81 We 
not only try a range of parameter values, but also check the 
effect of temporarily decoupling the hydrodynamical forces 
on wind particles. This section also investigates a promising 
way of injecting SN feedback in thermal form. Another form 
of feedback from young stars is studied in section IT9] where 
we discuss the results of various simplified implementations 
of radiatively driven winds. Finally, we study the effect of 
AGN feedback in section |4T0] 



4.1 Cosmology 

To investigate the dependence on cosmology and in order 
to facilitate comparisons to earlier work, we change the cos- 
molo gical parameters fro m our fiducial WMAP 3-year val- 
ues (jSpergel et all |2007| ) to the cosmology used in many 
studies including the Millennium Simulation (|Springel et al.l 
12005b). We refer to this latter set of cosmological parame- 
ters, which were chosen to be consistent with a combined 
analysis of the 2-degree field galaxy redshift survey and the 
first-year WMAP data, as the 'Millennium cosmology' and 
denote the models MILL. The Millennium cosmology uses 
the cosmological parameter values {fim, fib, JIa, erg, ris, /i} = 
{0.25,0.045,0.75,0.9,1.0,0.73}. Note that because of the 
change in the values of Sim and ilb, the dark matter and 
the (initial) baryonic particle masses are, respectively, 4.5 
and 7.7 percent higher for the MILL run than for the REF 
model. 
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Table 2. List of main physics variations employed in the OWLS project. From left to right the columns show the simulation name, 
whether or not the simulation was run in the 25 Mpc//i and 100 Mpc//i boxes respectively, the section number containing the description 
of the model, and a very brief description of the changes in the model relative to the REF simulation. Except for the MILL runs, all 
simulations that were run using the same box size used identical initial conditions. 



Simulation 


L025 


LlOO 


Section 


Description 


AGN 


V 


V 


14 101 


Includes AGN 


DBLIMFCONTSFV1618 


V 


V 


14.7.21 


Xop-heavy IMF at high pressure, cont. SF law, extra SN energy in wind velocity 


DBLIMFV1618 


V 


V 


14.7.21 


Xop-heavy IMF at high pressure, extra SN energy in wind velocity 


DBLIMFCONTSFMLI4 


V 


V 


14.7.21 


Xop-heavy IMF at high pressure, cont. SF law, extra SN energy in mass loading 


DBLIMFMLI4 


V 


V 


14.7.21 


Xop-heavy IMF at high pressure, extra SN energy in mass loading 


EOSlpO 


V 


V 


l44l 


Slope of the effective EOS changed to ^eff ~ f 


EOSlp67 


V 






Slope of the effective EOS changed to 'ycff — 5/3 


IMPS ALP 


V 


J 

V 


14.7.11 


Salpeter (1955) IMF 


IMFSALPMLl 


V 




|4.7.1| 


Salpeter (1955) IMF; wind mass loading Vj = 2/1.65 


MILL 


V 


V 




Millennium simulation cosmology, Tj ^ A (twice the SN energy of REF^ 


NOAGB.NOSNIa 




J 

V 


11:61 


No mass loss from AGB stars and SNIa 


NOHeHEAT 


J 

V 






No extra heat input around helium reionization 


NOREION 


V 




[431 


No hydrogen reionization 


NOSN 


J 

V 


J 

V 


01 


No SN energy feedback from SNe 


NOSN.NOZCOOL 


J 

V 


J 

V 


11:21 


No SN energy feedback from SNe and cooling assumes primordial abundances 


NOZGOOL 


J 

V 


J 

V 


112] 


Cooling assumes primordial abundances 


REF 


J 

V 


J 

V 


El 


Reference model 


REIONZ06 


V 




11731 


Hydrogen reionization occurs at z = Q 


REIONZ12 


J 

V 




1131 


Hydrogen reionization occurs at 2 = 12 


SFAMPLx3 




- 


|4.5.2| 


Normalization of Kennicutt-Schmidt SF law increased by a factor of 3 


SFAMPLx6 






|4.5.2| 


Normalization of Kennicutt-Schmidt SF law increased by a factor of 6 


SFSLOPElp75 






14.5.21 


Slope of Kennicutt-Schmidt SF law increased to 1.75 


SFTHRESZ 






|4.5.1| 


Critical density for onset of SF is a function of mctallicity (Eq. |1J 


SNIaGA USS 






HD 


Gaussian SNIa delay function 


WDENS 






|4.8.1| 


Wind mass loading and velocity depend on gas density (SN energy as REF) 


WHYDRODEC 






|4.8.2| 


Wind particles are temporarily hydrodynamically decoupled 


WML1V848 






|4.8.1| 


Wind mass loading 77 = 1, velocity t)„ = 848km/s (SN energy as REF) 


WML4 






11:81 


Wind mass loading 77 = 4 (twice the SN energy of REF) 


WML4 V424 






|4.8.1| 


Wind mass loading rj = 4; wind velocity = 424km/s (SN energy as REF) 


WML8V300 






|4.8.1| 


Wind mass loading 77 = 8; wind velocity -Uw = 300km/s (SN energy as REF) 


WPOT 






11:91 


Wind mass loading and vel. vary with grav. potential ("Momentum-driven") 


WPOTNOKIGK 






11:91 


Same as WPOT except that no extra velocity kick is given to winds 


WTHERMAL 






|4.8.3| 


SN energy injected thermally 


WVCIRC 






11:91 


Wind mass loading and vel. vary with halo circ. vel. ("Momentum-driven") 



The main differences with respect to the reference 
model are the values of fib and ag which are, respectively, 8 
and 22 percent higher for MILL than for REF. Both changes 
are expected to increase the SFR. The higher value of as has 
a particularly large effect at high redshift, because structure 
formation proceeds faster in the MILL cosmology. In order 
to roughly match the peak in the observed SFH, we dou- 
bled the mass loading factor to ?7 = 4 for the SN driven 
winds. Hence, the winds account for 80 percent of the avail- 
able energy from SNe. To isolate the effect of cosmology, 
we therefore compare the MILL simulation to model WML4 
which employs the same wind parameters, but is otherwise 
identical to the reference model. 

Fig. [5] compares the SFHs in the MILL (dashed, red) 
and WML4 (dot-dashed, blue) runs. The change from the 
WMAP-3 to the MILL cosmology strongly boosts the SFR. 
The difference increases with redshift from about 0.2 dex 
at z = to 0.34 dex at 2; = 2 (for both box sizes). By 
z — 9 the difference has increased to 1.0 dex. Clearly, for a 
quantitative comparison with observations, it is important 
to use the correct cosmology. At high redshift, when the 



haloes that dominate the SF in the simulation correspond 
to rare fluctuations, the predicted cosmic SFR becomes ex- 
tremely sensitive to the value of erg. We will show in Haas 
et al. (in preparation) that the differences are much smaller 
for haloes of a fixed mass, which implies that the change in 
the halo mass function accounts for most of the differences 
in the SFHs predicted for the two cosmologies. 

The olive, dotted curve in Fig. [S] shows the SFH pre- 
dicted by the Gala xies-Intergalactic Medium Interaction 
Calculation (GIMIC. [Crain et al.|[2009l ). GIMIC consists of 
a series of hydrodynamical simulations that zoom in on 
25 Mpc subvolumes of the 500 Mpc dark matter only 
Millennium simulation and was run using the same code 
and parameter values as MILL. Fig. [S] shows the SFH com- 
puted from the weighted average of the five GIMIC sub- 
volumes. The particle mass (gravitational softening) used 
for the GIMIC runfl is 8 (2) times larger than that of our 
25 Mpc box and thus 8 (2) times smaller than for our 

^ Xhese are the numbers for the intermediate resolution GIMIC 
runs. Xhe high-resolution runs use the same particle masses and 
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Figure 5. The effect of cosmology on the cosmic SFH. The curves 
show the cosmic SFR density as a function of redshift (bottom 
X-axis) and lookback time (top x-axis) for models MILL (red, 
dashed) and WML4 (blue, dot-dashed). For comparison, the re- 
sults for the reference model are also shown (black, solid). Results 
are shown for both the 25 and the 100 Mpc boxes, with the 
smaller box predicting higher SFRs at high redshift. All simula- 
tions used 2 X 512^ part icles. The data points show the compila- 
tion of observations from lHopkins fc Beacornl | |2006|) , converted to 
our IMF and our fiducial cosmol ogy. For compariso n, we also show 
the SFH predicted by GIMIC l lCrain et aLlboOOl. olive, dotted). 
Models MILL and WML4 differ only in terms of the assumed 
cosmology. In particular, model MILL assumes values for Ht, and 
as that are, respectively, 8 and 22 per cent higher than our fidu- 
cial values. Both MILL and WML4 assume that the wind mass 
loading, and hence the fraction of the SN energy that is injected 
in the form of winds, is twice as high as for model REF. Note 
that the lookback time axis (top x-axis) applies only to the REF 
cosmology. Although the data points assume the REF cosmology, 
they would be very similar for the MILL cosmology. Compared to 
our fiducial WMAP-3 cosmology, the MILL cosmology predicts 
higher SFRs, particularly at high redshift. 



100 Mpc box. At very high redshift the SFR in the 
GIMIC run is intermediate between that of our two MILL 
box sizes. For 2 < 7 it is very close to that of the 25 Mpc 
box and at 2 < 3 it agrees with the 100 Mpc run, al- 
though the GIMIC SFR falls of somewhat more steeply be- 
low 2 = 2. These differences are exactly what is expected 
from resolution effects as can be seen by comparing the SFHs 
for the 25, 50 and 100 Mpc boxes shown in the right 
panel of Fig. The excellent agreement confirms that our 
box sizes are sufficiently large to obtain a converged predic- 
tion for the cosmic SFH. 



4.2 Metal-line cooling 

Simulations without any radiative cooling are of interest for 
the study of hot gas in groups of clusters of galaxies (we have 
run such a simulation for this purpose in the 100 Mpc 
box), but in order to form stars, the gas must be able to 
radiate away its binding energy. Despite the importance of 
cooling, most cosmological studies still use highly simplified 



force resolution as our 25 /i ^ Mpc box, but they end at 
and do not include the highest density subvolume. 



Figure 6. As Fig.fS] but comparing the SFHs for models with and 
without metal-line cooling, both in the presence and absence of 
SN-driven winds. Except at very high redshift, metal-line cooling 
strongly increases the SFR. The boost due to metal-line cooling is 
greater when SN feedback is included, which implies that metals 
radiate away a significant fraction of the energy injected by SNe. 



prescriptions, ignoring metal-line cooling or including it un- 
der the assumption of coUisional ionization equilibrium and 
fixed relative abundances. Our simulations are the first to 
compute the cooling rates element-by-element and the first 
high-resolution simulations of cosmological volumes that in- 
clude the effect of photo-ionisation on the heavy elements. 

Fig-Elcompares the reference simulations with runs that 
ignored metal-line cooling {NOZCOOL; dashed, red). As ex- 
pected, the two agree at very high redshift where there has 
not been enough time to enrich the gas s ignificantly and 
where much of the gas falls i n cold (e.g. IWhite fc FrenkI 
1991!; iBirnboim fcP ckcl 2003; iKeres et a.1.1 12005| ). At late 
times, the runs without metal-line cooling consistently pre- 
dict lower SFRs. For our high-resolution L025N512 runs the 
difference increases with cosmic time to 0.3 dex at 2 = 2. 
While the SFR increases to 2 = 2 for REF, it peaks at 2 = 4 
when metal-line cooling is ignored. Interestingly, for the 
LlOO runs the difference decreases after peaking at about 
0.4 dex around 2 = 0.4. 

Also shown in Fig. [6] are foui[f| runs without SN- 
driven winds (but still including metal production and mass 
loss from SNe), both with and without metal- line cooling. 
Clearly, SN feedback strongly suppresses the SF, a point 
that we will come back to in section 

Interestingly, while metal-line cooling also enhances the 
SFR in the absence of SN feedback, its effect is smaller than 
when SN feedback is included. Put another way, the fac- 
tor by which SN feedback reduces the SFR is smaller when 
metal-line cooling is included. There are two possible ex- 
planations for this effect, which may both be right. First, 
metal-line cooling may reduce the efficiency of SN feedback, 
probably because it increases radiative losses in gas that has 
been shock-heated by the wind. Second, SN feedback may 
increase the effect of metal-line cooling, probably because it 
increases the fraction of the gas that is enriched. The former 
explanation is likely to be most relevant for galaxy groups. 



* Note that NOSN.NOZCOOL.L025N512 was stopped earlier. 



© 2007 RAS, MNRAS 000.[Tlf27l 



12 J. Schaye et al. 



as we will show elsewhere that galactic winds do not domi- 
nate the enric hment of the intrag r oup m edium. 

Recently, IChoi fc Nagamind (120091 1 have also investi- 
gated the effect of metal-line cooling on the SFH. While they 
also used GADGET, their simulations used about an order of 
magnitude fewer particles and were stopped at higher red- 
shifts. They did not include stellar evolution, they assumed 
a different cosmology and used different subgrid prescrip- 
tions for_SF_aBd^J£Jeedb^^ cooling rates were taken 
from [Sutherland fc Dopital (|l993l ) and assume fixed relative 
abundances and coUisional ionization equilibrium. Hence, 
there are many differences compared to our implementation 
of cooling, as we do include stellar evolution and compute 
the cooling rates element-by-element, including the effects 
of photo- ionisation by the meta-galactic UV/X-ray back- 
ground. For their simulations with resolutions similar to our 
L025N512 runs (but a much smaller box), they found that 
metal-line cooling increases the SFR by less than 0.1 dex at 
2 = 3, whereas we find 0.16 dex. In simulations with resolu- 
tion comparable to our L100N512, they found an increase of 
about 0.25 dex by 2 = 1 whereas we find 0.3 dex. Thus, the 
results are broadly consistent although we find a somewhat 
larger boost due to metal cooling. 

Clearly, except at very high redshift, metal-line cool- 
ing is very important. It boosts the SFRs by allowing more 
of the gas that accretes onto haloes to cool and by reduc- 
ing the efficiency of SN feedback. Without metal cooling, 
predictions for the total amount of stars formed could eas- 
ily be low by a factor of two. Note, however, that we may 
have overestimated the effect of metal-line cooling for mas- 
sive galaxies, because we have not included any feedback 
processes capable of stopping cooling fiows in such systems. 
This results in an overestimate of the SFRs and the metallic- 
ities in the central regions of groups and clusters of galaxies . 
On the other hand, as discussed in |w lersma et al.1 (l2009bl ). 
the fact that SFH underestimates small-scale metal-mixing 
(because metals are stuck on particles) causes us to under- 
estimate the total mass that has been enriched, while over- 
estimat ing the metallic i ty of th e particles that have received 
metals. IWiersma et al.l l|2009bD found that the net effect of 
increased metal mixing is to boost the SFR. 



4.3 Reionization 

As our simulations do not include radiative transfer, we 
need to assume the background radiation is uniform and 
that the gas is optically thin. Ifydrogen reionization is 
thus implemented in ou r simulations by switching on the 
iHaardt fc Madaul (|200ll ') model for the ionizing background 
radiation at some redshift Zr, corresponding to the epoch of 
reionization. Note, however, that we assume that a photo- 
dissociating background is already present ai z > Zr, which 
effectively s uppresses molecular coo ling at all redshifts. As 
described in |w lersma et al.1 (|2009al ). switching on the ion- 
izing radiation results in a sudden increase in the radiative 
heating rate and a sudden decrease in the radiative cooling 
rate above 10^ K. As a result, cold gas is quickly heated to 
T ~ 10'* K, removing gas from haloes with virial temp era- 
tures < 10^ K (e.g. Couchman fc Rccs 1986; Okamoto o t al.l 
l2008l ). 

Fig.Qcompares L025N512 runs with Zr = 12, 9 (i.e. the 
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Figure 7. As Fig. [5] but comparing the SFHs for models in which 
hydrogen was reionized — by switching on a uniform ionizing 
background — at redshift = 12 {REIONZ12), 9 (REF), 6 
{REIONZ06), or not at all (NOREION). All simulations use a 
25 Mpc box and 2 X 512^ particles. After the ionizing back- 
ground is switched on, the SFR quickly changes to the rate pre- 
dicted by the model with the highest redshift of reionization. Note 
that the factor by which reionization suppresses the SFR is lim- 
ited by the resolution of the simulations and will be more severely 
underestimated at higher redshifts. Turning off the extra heat in- 
put of 2 eV per atom around helium reionization (2 3.5, model 
NOHeHEAT) has no discernible effect on the SFH. 



reference model), 6 and (i.e. no reionization After reion- 
ization, the SFR deviates from the curve corresponding to 
the simulation without an ionizing background (NOREION) 
and quickly asymptotes to the model with the highest red- 
shift of reionization {REI0NZ12). Apparently, the gas in the 
simulation rapidly l oses memory o f the t ime of reheating, as 
was also found bv iPawUk et al.l (|2009l '). This is expected, 
as the sound-crossing time scale is only 10*yr(//l kpc) for 
10" K gas. 

Naively, one would have expected the suppression of 
the SFR due to photo-ionization to decrease at late times, 
as haloes with TVir ^ lO** K start to dominate the cosmic 
SFR. Interestingly, we do not find this. If anything, the sup- 
pression keeps increasing with time, reaching 0.15 dex (a 
thirty percent reduction) by z = 2. While this is proba- 
bly mostly a resolution effect, it does indicate that photo- 
ionization also reduces the SFR in haloes with higher virial 
tempe ratures, either because of the reduction of the cooling 
rates (jWiersma et al.ll2009ah or because i t makes the cold 
gas m ore susceptible to galactic winds l|Pawlik fc Schav3 
1 200^ . 

We emphasize that because of our limited resolution 
(Haas et al., in preparation, show that we underestimate the 
SFR in haloes with masses less than 10^" Mq) and because 
we assume the presence of a photo-dissociating background 
at all redshifts, it is likely that we have strongly underesti- 
mated the reduction of the SFR due to photo-heating and 
that this underestimate becomes more severe at higher red- 
shifts. Moreover, our assumption that the gas is optically 



3 The run with Zr = 12 uses the z = 9 iHaardt fc Madaul 1I2OOII') 
model for 2 = 9— 12. 
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Figure 8. As Fig. [5] but comparing tlie SFHs for models assum- 
ing different equations of state for the unresolved ISM. The slope 
of the polytropic EOS imposed on the ISM is 1.0 (i.e. isothermal) 
for model EOSlpO (red, dashed), 4/3 for REF (solid, black), and 
5/3 (i.e. adiabatic) for EOSlp67 (blue, dot-dashed). The SFH is 
insensitive to the assumed slope of the polytropic EOS that is 
imposed onto the ISM. 



thin results in an underestimate of the heating rates during 
reionization. 

Helium is thought to have been reionized around z ~ 3.5 
and the increase in the photo-heating rates associated with 
this event can explain the relatively high temperature of the 
IGM inf erred from observa tions of quasar a bsorption spec- 
tra (e.g . ISchave et ahlbOOOl ) . As described in |w iersma et al.l 
lj2009bl ). by injecting 2 eV per atom at « « 3.5, we are able 
to match the observationally inferred temperatures. Fig. [7] 
shows that omitting this extra heat (NOHeHEAT) does not 
yield any noticeable changes in the SFH. This is not surpris- 
ing, as the temperature increase is confined to low-density 
gas, far away from galaxies, for which adiabatic cooling dom- 
inates over radiative cooling. 



4.4 The equation of state of the ISM 

Our simulations have neither the resolution nor the physics 
to model the multiphase ISM. As discussed in section 13.1.31 
we therefore impose a polytropic EOS with slope 7etF = 4/3 
for gas with densities that exceed our SF threshold of 
riH = 0.1 cm~'^. This slope was chosen because it results 
in a constant Jeans mass and thus suppresses artificial frag- 
mentation. In this se ction we will check the effect of varying 
the slope of the EOS. ISchave fc Dalla Vecchial (j2008h showed 
that changes in the EOS can significantly alter the morphol- 
ogy of galaxies. A softer EOS results in tighter spiral arms, 
thinner disks, and increased fragmentation. 

Different groups use different pr e scriptio ns for the ISM. 
For example, ISpringel fc Hernguistl (|2003al ) use a compli- 
cated function that results from a semi-analytic model of 
the multiphase ISM. They interpret the pressure implied by 
their EOS, which is steeper than 4/3 at densities similar to 
our SF threshold, as a form of SN feedback. In the past, 
many cosmological simulations have been run that do not 
impose an EOS, but which also do not include the physics 
necessary to model the cold interstellar gas phase (e.g. ra- 



Figure 9. As Fig. \5\ but comparing the SFHs for the refer- 
ence model (black, solid), which uses a fixed threshold density 
for SF, to that of model SFTHRESZ (red, dashed), for w hich the 
SF threshold decreases with metallicity as predicted by ISchavd 
(2004). The threshold densities in the two models agree for a 
metallicity of 0.1 Zq. Both simulations use a 25 Mpc box 
and 2 X 512^ particles. At very high redshift the metallicity is low 
and the total SFR is smaller for SFTHRESZ because it has a 
higher threshold density at this point. Below z = 6 the situation 
is reversed, but the difference between the SFHs is very small, 
suggesting that the cosmic SFR is dominated by galaxies that 
are able to regulate their SFRs. 



diative transfer and molecule formation). Such simulations 
effectively use an isothermal EOS. 

Fig. |8] compares the SFHs of runs with j^ff = 1 (i-e. 
isothermal), 4/3 (REF), and 5/3 (i.e. adiabatic). Clearly, 
changes in the EOS do not have a significant effect on the 
predicted SFH. This may be surprising, given that the EOS 
can strongly affect the structure of galaxies. 

One reason why the results are insensitive to 
the EOS is that we use t he prescription for SF of 
ISchave fc Dalla Vecchial (j2008l ). As discussed by these au- 
thors, for self-gravitating systems such as galaxies, the ob- 
served Kennicutt-Schmidt surface density law is in effect a 
pressure law. By implementing it as a pressure law, we can 
thus reproduce the observed SF law independently of the 
assumed EOS of the star-forming gas. Previous cosmologi- 
cal simulations have, however, used volume density laws, in 
which case the predicted Kennicutt-Schmidt law must de- 
pend on the assumed EOS because the latter sets the scale 
height of the disk. If the EOS is changed, then the same sur- 
face density corresponds to a different volume density, but 
the relation between surface density and pressure will remain 
unchanged. It is therefore not clear whether the results of 
previous simulations are as insensitive to the imposed EOS 
as we find here. However, we will show in the next section 
that the results are in fact insensitive to the gas consumption 
time scale, because SN feedback enables galaxies to regulate 
their SFRs. 



4.5 Star formation 

4-5.1 The star formation threshold 

ISchavd l|2004 ) argued that there is a critical density for the 
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formation of a cold, interstellar gas phase and that the tran- 
sition from the warm to the cold gas phase triggers gravita- 
tional instability on a wide range of length scales. Gas with 
densities below the threshold is kept warm (T ~ 10^ K) and 
stable by the presence of a UV background. The predicted 
critical gas surface density Eg ~ 3 — 10 Mqpc"^, which 
agrees well with SF thresholds inferred from Ha observations 
of nearby galaxies, corresponds to nn ~ 10~^ — 10~^ cm~'^ 
for a self-gravitating dis k at 10"^ K. Ou r reference model uses 
riH = 10"^cm"3. The ISchavj (|2004l) model predicts that 
the critical density for SF is a weakly decreasing function of 
metallicity. We have therefore run a simulation, SFTHRESZ, 
that uses th e predicted scaling (equations 19 and 24 of 
ISchavell200j . valid for Z = lO""* - 10 Zq), 
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n'^{Z) = 10"^ cm" 



Q.lZr, 



(4) 



where Z is the g as meta ll icity and we used Zq =0.02 for 
consistency with ISchavd (|2004h . If the metallicity is zero 
then we set — 10 cm~^. 

Fig. |9] compares models SFTHRESZ and REF. At very 
high redshift the metallicity is low and the threshold den- 
sity is higher than in the reference model. This results in 
a decrease in the SFR that drops rapidly from 0.3 dex at 
2 = 10 to zero by z; = 6. For z <G the SFR is slightly higher 
than in the REF model, which indicates that the metallic- 
ity of the star- forming gas is typically higher than 0.1 solar, 
but the effect is marginal. Apparently, after a brief period in 
which the SFR is dominated by haloes that are just resolved 
and therefore just starting to form stars, the predicted SFRs 
become insensitive to the SF threshold. This suggests that 
the galaxies are able to regulate their SFRs. We will provide 
more evidence for this below. 



4-5.2 The Kennicutt-Schmidt star formation law 



As discussed in Sec. 13.1.31 gas on the effective EOS is al- 
lowed to form stars at a pressure-dependent ra te that re- 
prod uces the observed Kennicutt-Schmidt law (|Kennicuttl 



1 19981 ). E. 



A(Sg/l Mopc"^)", with A = 1.515 x 



10"" M0yr"^kpc"^ and n = 1.4. The normalization (A) 
and slope (n) are constrained by ob servations, but remain 
controversial (e.g. iBlanc et al.ll2009l ). To develop an under- 
standing of the physical role of the SF law, we have carried 
out one run with a different slope and two with different 
amplitudes. 

Fig. [To] compares a run with n — 1.75 (model 5^5- 
LOPElp75; red, dashed) with our reference model, which 
uses n = 1.4. The SF laws are in both cases normalized at 
Eg = 1 Mqpc"^, which is below the threshold and hence 
implies that the SFR is higher for all densities in the run 
with the steeper slope SF law. For z > 6 the cosmic SFR 
is indeed higher in the run with n — 1.75. This is expected, 
because at these high redshifts the SFR is dominated by 
haloes that are just resolved and therefore just starting to 
form stars. These galaxies have not yet had time to become 
self-regulating and their SFRs are inversely proportional to 
the gas consumption time-scales implied by the SF law. 

Below z = 6, however, the SFRs in the two runs are 
nearly indistinguishable. This strongly suggests that the 
galaxies are regulating their SFRs such that they produce 
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Figure 10. As Fig. [S] but comparing the SFHs for models with 
varying Kennicutt-Schmidt SF laws. Model SFSLOPElplS (red, 
dashed) assumes a power-law slope n = 1.75 whereas the other 
models use our fiducial value n = 1.4. For models SFAMPLxS 
(blue, dot-dashed) and SFAMPLxd (olive, dotted) the amplitude 
of the SF law has been multiplied by factors of 3 and 6, respec- 
tively. All simulations use a 25 Mpc box and 2x512^ particles. 
At very high redshift, when the SFR in the simulations is domi- 
nated by poorly resolved haloes, a more efficient SF law yields a 
higher SFR. After this initial phase the SFH is insensitive to the 
assumed SF law, which suggests that it is dominated by galaxies 
that are able to regulate their SFRs. 



the same amount of stars, and thus the same amount of SN 
energy, irrespective of the gas consumption time scale. If a 
galaxy of a given halo mass, and hence with a fixed accre- 
tion rate, injects too little SN energy for a galactic outflow to 
balance the accretion rate, then the gas fraction, and hence 
the SFR, will increase. If, on the other hand, the SN rate is 
higher than required to balance the infall, then the gas frac- 
tion, and thus the SFR, will decrease. We thus expect that 
when the SF efficiency is changed, the galaxies will adjust 
their gas fractions so as to keep their SFR fixed. In Haas et 
al. (in preparation) we show that this is indeed what hap- 
pens. 

Finally, Fig. [10] shows that models in which the ampli- 
tude of the SF law is multiplied by factors of three [SFAM- 
PLxS; blue, dot-dashed) and six {SFAMPLxG; olive, dotted), 
respectively, show the same behavior. Initially, the SFR in- 
creases with A, but the SFR then quickly asymptotes to 
a fixed SFH. Observe that the two runs with higher am- 
plitudes converge to a common evolution before the refer- 
ence model joins them. This is because galaxies can reg- 
ulate their SF more quickly if the SF efficiency is higher. 
Apparently, the cosmic SFR in the reference model only be- 
comes dominated by self-regulated galaxies by z = 6. Note 
that higher resolution simulations may well find that self- 
regulation dominates already at higher redshifts because 
they can resolve SF in the progenitors of our lowest mass 
galaxies. 

4.6 Intermediate mass stars 

Previous numerical studies of the cosmic SFH have mostly 
used the instantaneous recycling approximatio n (but see e.g. 
lOppenheimer fc Dave|[2008l : ICrain et al.|[2009l ). which means 
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Figure 11. As Fig. [5] but comparing the SFHs for models that 
vary in their treatments of intermediate mass stars. In model 
NOAGB_NOSNIa (red, dashed) mass loss by intermediate mass 
stars and SN Type la has been turned off. While model SNIa- 
GAUSS (blue, dot-dashed) includes these processes, it assumes 
a Gaussian time-delay function for Type la SNe instead of the 
e-folding model used in the other simulations. All simulations use 
a 100 h~'^ Mpc box and 2 X 512^ particles. While the SNIa time 
delay function is unimportant, mass loss by AGB stars provides 
fresh fuel for SF and releases metals back into the ISM, thereby 
boosting the SFR at late times. 



that star particles eject all the products of stellar evolution 
immediately following their formation. Moreover, individ- 
ual elements are typically not tracked. Instead, each gas 
element carries only a single metallicity variable and rela- 
tive abundances are assumed to be solar. Furthermore, such 
simulations neglect mass loss, i.e., star particles change the 
metallicity of their neighbors, but not their masses. As dis- 
cussed in section 13.1.41 we follow the timed release of 11 
elements by intermediate mass stars (SNIa and AGB stars) 
and massive stars. 

To check the importance of intermediate mass stars, 
which eject much of their mass hundreds of millions to bil- 
lions of years after their formation and which are responsible 
for most of the mass lost by stellar populations, we have run 
two L100N512 simulations. In simulation NOAGB.NOSNIa 
we do not allow intermediate mass stars to release mass, 
leaving massive stars, which evolve on timescales of < lO'^ yr, 
as the only mechanism for releasing metals. To assess the im- 
pact of our choice of the distribution of SNIa progenitor life- 
times, we ran a simulation (SNIaGAUSS) that uses a Gaus- 
sian rather than an e-folding time delay funct ion. This was 
motiv ated by the high redshift observations of lPahlen et al.l 
l|2004l , which show a marked decline in the SNIa rate beyond 
z = I. The paramete rs of the delay model are a = 0.66 Gyr 
and r = 3.3 Gyr (see lWiersma et al.l[2009bl l . 

Fig. [TT] shows that the shape of the SNIa delay func- 
tion does not have a significant effect on the predicted SFIf. 
Turning off both mass loss by SNIa and AGB stars results, 
however, in a strong reduction of the SFR at late times. 
While the reduction factor is still very small at z = 2, it in- 
creases steadily thereafter to about 0.21 dex at z = 0, which 
corresponds to a 40 percent reduction. Given that the SNIa 
delay function does not matter, the difference must come 
mostly from mass loss by AGB stars. It is not important 



Figure 12. As Fig. [5] but comparing the SFHs for models as- 
suming a Salpeter IMF (IMFSALP; red, dashed) to the reference 
model (black, solid), which assumes a Chabrier IMF. The ampli- 
tude of the SF law is taken from observations and has therefore 
been rescaled to the assumed IMF. It is a factor of 1.65 higher for 
the Salpeter IMF. Model IMFSALPMLl assumes a Salpeter IMF 
and uses a wind mass loading factor that is a factor of 1.65 smaller 
(i.e. r) = 2/1.65) than that used in the other models, which ac- 
counts for the change in the number of SNe per unit stellar mass. 
Note that the observed data points assume a Chabrier IMF. They 
need to be shifted upwards by a factor 1.65 (0.22 dex) to compare 
with models assuming a Salpeter IMF. Initially the SFR scales 
with the amplitude of the SF law and a Salpeter IMF produces a 
higher SFR. Later on the SFR is smaller for a Salpeter IMF be- 
cause of the decreased importance of metal-line cooling (because 
less metals are produced and a greater fraction are locked up in 
stars) and stellar mass loss. However, the smaller number of SNe 
per unit stellar mass more than compensates for this effect, at 
least for z > 2. 



before z = 2 because there has not been sufficient time for 
a substantial fraction of the stars to reach the AGB phase. 
Note, however, that higher resolution simulations will pre- 
dict higher SFRs at high redshift and may therefore find 
that AGB mass loss becomes important earlier. 

Mass loss by AGB stars provides fresh fuel for SF and 
releases metals that were locked up in stars. This reduces the 
sharpness of the drop in the SFR with time, worsening the 
agreement with observations. Simulations that ignore this 
process, will overestimate the steepness of the drop following 
the peak in the cosmic SFR. 



4.7 The stellar initial mass function 

4.7.1 A Salpeter IMF 

Our reference model assumes a Chabrier IMF, but much of 
the literature uses a Salpeter IMF. The two IMFs have simi- 
lar shapes above 1 Mq , but while the Salpeter IMF is a pure 
power-law, the Chabrier IMF includes a lognormal decrease 
at the low mass end which results in a much lower stellar 
mass-to-light ratio. Because most of the ejected metal mass 
and all of the energy from core collapse SNe is produced by 
massive stars, the Salpeter IMF is less efficient in enriching 
the gas and driving outfiows per unit stellar mass formed. 

In order to assess the effect of changing the IMF we 
ran a simulation employing a Salpeter IMF {IMFSALP), 
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using the same range of stellar masses as we used in the 
reference model (i.e. 0.1 — 100 M©). This simulation used 
Kennicutt's original normalization for the amplitude of th e 
SF law {A = 2.5 x 10"* Mr;, vr~^ kpc~^: lKennicut^ll998l l. 
as he assumed the same IMF. Recall that this value is a 
factor 1.65 greater than the amplitude assumed in REF (see 
Section I3.1.3P . 

Fig. HH compares the SFHs predicted by the two IMFs. 
Initially the Salpeter IMF gives a slightly higher SFR be- 
cause of the higher SF efficiency implied by the change in the 
SF law (see also section |4. 5. 2p . However, after a short initial 
phase the SFR falls below that of the reference model. By 
z = the difference has increased to 0.2 dex for L100N512. 
This behavior can be explained by the fact that a Salpeter 
IMF produces less metals and returns less mass (and thus 
releases less metals that were locked up in stars) per unit 
stellar mass formed. Hence, metal-line cooling is less effi- 
cient for a Salpeter IMF. Indeed, the predicted SFRs fall in 
between those for the reference model and the run without 
metal-line cooling (c.f. Fig. |B} . 

However, a Salpeter IMF does not only produce less 
metal mass per unit stellar mass, it also produces fewer 
SNe. Assuming that the total energy in SNe scales as the 
total number of ionizing photons, the difference is a factor 
of 1.65 (Section EUnj. Thus, model IMFS ALP uses 66 per- 
cent of the SN energy to drive winds, whereas REF used 
only 40 percent. For consistency, we therefore ran another 
L025N512 simulation, model IMFSALPMLl, that is identi- 
cal to IMFSALP, except that the wind mass loading factor 
was reduced by a factor 1.65 to 77 = 1.2. Fig. 1121 shows that, 
as expected, this run yields a higher SFR than IMFSALP, 
although the two converge for z > 9 where there has not 
been sufficient time for the simulated galaxies to regulate 
their SF. In fact, with this change, a Salpeter IMF yields a 
higher SFR than a Chabrier IMF. 



4-1-2 A top-heavy IMF at high pressures 

Observational determinations of the IMF are extremely dif- 
ficult. In particular, extragalactic observations are usually 
only sensitive to the light emitted by massive stars, either 
directly or indirectly via dust grains. While the IMF is usu- 
ally assumed to be universal, it is expected to be top-heavy 
(or bottom-hght) at very high redshift and low metallicity 
(e.g. iLarsonlFiggsl ) and both observations and theory suggest 



tic center and starburst galaxies (e.g. 


Padoan et al. 


1997 


BauKh et al.l 20051: Klessen et al. 2007; 


Maness et al. 


2007 


DabrinEhausen et al.l 20091: Bartko et al.l 2009t). 



We have performed a series of runs to investigate the 
possible effects of an IMF that is top-heavy at high pres- 
sures, as may be the case in starbursts and in galactic 
nuclei. For simplicity, we assume the IMF switches sud- 
de nly from Cha b rier t o the top-heavy power-law proposed 
bv iBaugh et al.l l|2005l '). dN/dM oc (as compared to 

oc M~^-^ for the high- mass tail of the Chabrier IMF). 
The transition is assumed to take place at the pressure 
P/k — 2.0 X 10® cm~^ K (evaluated at the resolution limit 
of the simulations) , which was chosen because for this value 
~ 10~^ of the stellar mass in our simulations forms at higher 
pressures. This ensures that the top-heavy IMF is impor- 
tant, but not dominant. Of course, a discontinuous depen- 



dence on pressure is not physical, but it is simple and serves 
to illustrate the qualitative effects of a top-heavy IMF in 
starbursts. 

Assuming that the SN energy scales with the emissivity 
in ionising photons, our top-heavy IMF yields 7.3 times more 
SN energy per unit stellar mass formed. We have therefore 
increased the energy injected into the galactic wind by the 
same factor for star particles born out of high-pressure gas. 
In terms of our kinetic prescription for winds, we can either 
increase the mass loading or the wind velocity. We have tried 
both. Models MLI4 use a 7.3 times larger mass loading, 
while models VI 61 8 use a -/T^ times higher wind velocity. 

If the actual SF law were continuous with pressure, then 
a sudden change in the IMF would imply a sudden change 
in the rate of formation of massive stars, which would man- 
ifest itself as a discontinuity in the apparent SF law inferred 
from observations under the assumption of a universal IMF. 
However, the obse rved SF law appears t o be a continuous 
power- law (though iKrumholz et al.ll2009l suggest that there 
may be a kink at Eg ~ IO^Mqpc"^, which corr esponds 
rough ly to the pressure (see equations 20 and 21 of ISchayg 
l2004 l at which we switch IMFs). We therefore tried two 
possibilities: models DBLIMF (DBLIMFCONTSF) assume 
a continuous (discontinuous) rate of formation of massive 
stars, but a discontinuous (continuous) SF law. 

The left panel of Fig. [13] compares the SFHs of 
all four models for each of the two box sizes. Note 
that models DBLIMFML14JL100N512 and DBLIMFCON- 
TSFML14_L025N512 were stopped earlier than the other 
runs. The SFHs agree at early times, when SF is confined to 
low-mass haloes for which the gas pressure remains low. The 
models in which the extra SN energy was used to increase 
the wind mass loading predict SFHs that are similar to that 
of the reference model. This could mean that increasing the 
wind mass loading does not strongly boost the efficiency of 
SN feedb ack at high gas pressures . This agrees well with the 
results of D al ia Vecchia fc Schavel(|2008l ). who found that at 
high pressures the kinetic feedback becomes inefficient due 
to gas drag and that the pressure above which this occurs 
increases with the wind velocity. Indeed, as can be seen from 
Fig. 1131 the models in which the wind velocity is increased 
for the top-heavy IMF do show a strong reduction in the 
SFR. 

However, a top-heavy IMF not only yields more SN en- 
ergy, but also more metal mass per unit stellar mass formed. 
The associated increase in the metal-line cooling rates will 
boost the SFRs (see Section |4]2]). The relatively small differ- 
ence between the MLI4 and REF model could therefore also 
mean that the increased wind mass loading compensates for 
the higher cooling rates. 

The differences between DBLIMF and the correspond- 
ing DBLIMFCONTSF runs are small. This was to be ex- 
pected, as we already showed in section l4.5.2l that the SFHs 
are insensitive to the SF law because galaxies regulate their 
gas fractions so as to produce the same amount of SN energy, 
irrespective of the assumed SF efficiency. 

To compare with the observed data points, which were 
derived from observations of massive stars under the as- 
sumption of a Chabrier IMF, we have to multiply the rate 
of SF in the top-heavy mode by a factor 7.3. The right 
panel of Fig. [13] shows that doing so reduces the drop at 
late times, i.e. the SFR inferred under the assumption of 
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Figure 13. As Fig. [5] but comparing the SFHs for models with a top-heavy IMF in starbursts. The transition to a top-heavy IMF 
{dN/dM Qc compared to oc M~^'^ for Chabrier) happens suddenly at a pressure P/k = 2 X 10*^ cm-3 K. The top-heavy IMF 

produces 7.3 times more core collapse SNe per unit stellar mass formed. Models ML14 and V1618 use this extra energy to increase the 
wind mass loading and velocity, respectively. Models DBLIMFCONTSF assume a continuous SF law, whereas models DBLIMF assume 
the rate of formation of massive stars to be continuous. The left panel shows actual SFRs, whereas the SFRs have been rescaled to 
the ones that would be inferred under the assumption of a Chabrier IMF in the right panel. Comparisons to the observed data points, 
which assumed a Chabrier IMF, are only self-consistent for the right panel. The models start to differ when some galaxies have become 
sufficiently massive to form a fraction of their stars with a top-heavy IMF. Models with a top-heavy IMF form less stars, which indicates 
that the relative increase in the SN rate is more important than the increase in the metal production rate. Whether the SF law is 
continuous or not is unimportant. Using the extra SN energy to increase the wind mass loading is less effective than increasing the wind 
velocity. A top-heavy IMF in starbursts reduces the SFR in massive galaxies, but the effect on the formation rate of massive stars is 
much less strong. 



a universal IMF falls off less steeply than the actual SFR. 
Observe that the differences between models DBLIMF and 
DBLIMFCONTSF are also reduced, particularly for z < 1. 
This supports our proposal that because galaxies regulate 
their SF, they inject a fixed amount of SN energy for a given 
halo mass. 

Finally, we note that the agreement between the 25 and 
100 h-^Mpc boxes is much poorer for the models with a 
top-heavy IMF in starbursts than it was for the other mod- 
els. This reflects our choice to make the IMF a function of 
the pressure, which is 1-1 related to the gas density in our 
simulations since star-forming gas follows a polytropic EOS. 
Increasing the resolution decreases the mass above which 
haloes contain enough particles to sample the high-density 
tail of the gas distribution. Hence, lower halo masses will be 
able to form some fraction of their stars with a top-heavy 
IMF and thus suppress subsequent SF. Clearly, using pre- 
scriptions for feedback that are functions of density or pres- 
sure will make the results more prone to resolution effects. 

We conclude that although our toy models are too sim- 
ple and sensitive to resolution, it is clear that a top-heavy 
IMF in starbursts can serve to suppress SF in high mass 
haloes. This can result in a steeper drop in the SFR at late 
times, as suggested by observations (although the effect is 
less strong when only the rate of formation of massive stars 
is considered). The stronger suppression in massive haloes 
can also shift the peak in the SFH to higher redshifts. 

4.8 SN-driven winds 

It is well known that simulations without galactic winds suf- 
fer from a severe overcooling problem: the SFR greatly ex- 
ceeds the observational constraints and is usually only lim- 



ited by numerical resolution (e.g. iBalogh et al.|[200ll ). As 
Fig. [S] shows, except at very high redshift when the SFR 
in the real Universe is expected to be dominated by haloes 
that are below the resolution limit of our simulations, the 
runs without SN feedback do indeed produce far too many 
stars. Moreover, comparison of the SFHs for models REF 
and WML4 in Fig.[5]shows that doubling the SN energy that 
is injected, in this case by doubling the wind mass loading, 
further reduces the SFR. 

In sections 14.5.21 and |4T21 we showed that the SFH is 
insensitive to the assumed SF law. We concluded from this 
that SF in galaxies is self-regulating: the SFR adjusts such 
that outflows driven by feedback from massive stars balance 
the infall driven by gas accretion onto haloes and radiative 
cooUng. If the SF law is changed, then galaxies simply adjust 
their gas fractions in order to inject the same amount of SN 
energy into haloes of a given mass. If galaxies do indeed reg- 
ulate their SFRs in this manner, then we would expect the 
rate of energy injection into the winds to remain constant 
if the fraction of the SN energy that is injected is varied. 
In other words, the SFR should be inversely proportional 
to this fraction. We can test this by comparing model REF 
to simulation WML4, which injects twice as much SN en- 
ergy per unit stellar mass. Fig. [S] shows that while the SFR 
is indeed lower for WML4, the difference is always smaller 
than 0.25 dex, whereas we would have expected 0.30 dex 
(i.e. a factor of two) at late times, when the SFR is dom- 
inated by galaxies that have had enough time to become 
self-regulating. 

There are, however, good reasons why we would not 
expect the SFR to be exactly inversely proportional to the 
efficiency of the SN feedback. First, a change in the SFR 
does not only change the rate of energy injection into winds. 
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it also changes the rate of metal injection and the rate at 
which mass loss from AGB stars supplies the galaxy with 
fresh fuel for SF. Both of these effects would, however, tend 
to lower the amount of SN energy that is needed for self- 
regulation, which means that we would expect the SFR to 
decrease faster than linear with the SN efficiency. This is 
opposite to what is actually happening. Second, a reduction 
in the SFR implies a reduction in the stellar mass and hence 
in the gravitational force. Again, this would imply that less 
SN energy is required for self-regulation, which would lead 
to an even larger drop in the SFR, contrary to what the 
simulations predict. 

The reason why the SFR varies more slowly with the SN 
efficiency than the inverse proportionality we would naively 
expect, is likely that SN winds are not effective in high mass 
galaxies, at least when kinetic feedback is used with a ve- 
locity of 600 kms~^. If the feedback is inefficient, then we 
would not expect the galaxies to be able to self-regulate. In- 
deed, we will show in Haas et al. (in preparation) that, at 
2 = 2, the SFR in WML4 is in fact half that of REF for 
haloes with total mass less than 10^^ M0 and that the SN 
feedback becomes inefficient for higher halo masses. 

4.8.1 Varying the parameters at constant wind energy 

As discussed in section 13.1.51 we inject the energy from 
SNe in kinetic form. Newly formed star particles kick their 
gaseous neighbours with a constant velocity in a ran- 
dom direction. On average, the mass kicked is r] times the 
mass of the star particle. While the product r)v1, determines 
the energy of the winds and is therefore constrained by the 
energy available from SNe, it is not clear a priori what val- 
ues should be chosen for the individual parameters. Note 
that they cannot be taken directly from observations be- 
cause the parameter values refer to the properties of the 
wind at the inter-particle distance (neighbours of new stars) 
which varies and will typically not agree with the scales rele- 
vant for the observations. The observational constraints are 
usually inferred from the velocity offset and column den- 
sities of blueshifted absorption lines, but it is unclear at 
what distance from the source the absorption occurs (e.g. 
IVeilleux et ahlbOOSl ). Moreover, the absorption lines probe 
only the cold part of the outflow. The constraints on the ve- 
locity and mass loading of the hot wind are also very poor. 

Given the lack of observational constraints, one would 
hope that the results are insensitive to the amount of 
mass that a fixed amount of S N energy is distributed over. 
iDalla Vecchia fc Schavd l|2008h showed that this is indeed 
the case for the SFRs provided the galaxies are well resolved 
and the wind velocity exceeds a critical value that increases 
with the pressure of the ISM and thus also with halo mass. 
If, however, the wind velocity is too low, then the wind parti- 
cles are immediately stopped by drag forces and never leave 
the ISM. If the disks are unresolved, then the hydrodynamic 
drag is underestimated and for sufficiently low resolutions all 
particles that are kicked are able to escape the ISM. 

The runs with a top-heavy IMF in starbursts (Fig. I13p 
show that SF is much more efficiently suppressed if the extra 
SN energy (relative to a Chabrier IMF) is used to increase 
the wind velocity Hw than if it is used to increase the mass 
loading factor 77. Since in these models the feedback energy 
is only boosted in high-pressure gas, this suggests that the 
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Figure 14. As Fig. [5] but comparing the SFHs for models that 
all inject the same amount of SN energy per unit stellar mass (i.e. 
rjv'^ is constant), but that assume different wind velocities. Mod- 
els WMLxVyyy assume a wind mass loading r] = x and velocity 
''w = yyy kms~^. Model WDENS assumes that the wind veloc- 
ity scales with the local sound speed, which implies fw oc p^^^ for 
our EOS, normalized to the value used in the reference model at 
the SF threshold. At high redshift, when the SFR is dominated 
by poorly resolved, low-mass haloes, the SN feedback is more effi- 
cient if it is distributed over more mass. However, at low redshift 
the situation is reversed. This can be explained if the feedback 
becomes inefficient when the velocity falls below a critical value 
that increases with galaxy mass. Scaling the wind parameters 
with local properties, such as the density for model WDENS, can 
help to keep the feedback efficient, but it also makes the results 
sensitive to the resolution. 



wind velocity of 600 kms~^ that was used in the reference 
model is insufficient at the high pressures that we required 
for the IMF to become top-heavy. This implies that our 
default prescription for SN feedback is inefficient in high 
mass galaxies, which could account for the fact that the 
SFR drops off less rapidly at late times than is observed. 

To further investigate the dependence on the two indi- 
vidual wind parameters, we have run a series of simulations 
that all inject the same amount of SN energy per unit stellar 
mass as the reference model, but assume mass loading fac- 
tors that differ by factors of two, ranging from 1 to 8 in the 
25 /i"^ Mpc box {v„ oc 77"^/^ varies from 848 to 300 kms"^) 
and from 1 to 2 in the 100 Mpc box («„ varies from 848 
to 600 kms"^). Fig. [14] compares the SFHs of these runs. 
Clearly, the results are not just determined by the total en- 
ergy injected into the wind, which is identical for all the 
runs. At high redshifts the SFHs are similar, although the 
feedback is slightly more efficient for higher values of rj. How- 
ever, at late times the different SFHs start to diverge, with 
higher wind velocities suppressing the SF more strongly. 

These results are consiste nt with the conclusions of 
iDaUa Vecchia fc Schavd (|2008l ). As the universe evolves, 
stars typically form in more massive haloes and the mini- 
mum wind velocity for which the feedback remains effective 
thus increases. Hence, the redshift for which the feedback 
becomes inefficient decreases with increasing wind velocity. 
At early times, when many stars form in haloes that are 
poorly resolved, higher mass loading factors are more effi- 
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cient because all particles that are kicked from poorly re- 
solved haloes are able to escape the ISM. 

The SFH, including the redshift at which it peaks, is 
clearly sensitive to the poorly constrained parameters -q and 

The same is likely to be true for other types of subgrid 
prescriptions than kinetic feedback. Thus, unless one varies 
the parameters of the wind model, which is unfortunately 
not always done in the literature, one risks overinterpreting 
the results. 

If the main goal were to reproduce the observed SFH 
and if one were willing to accept the lack of "ab initio pre- 
dictive power" displayed by Fig. 1141 then one could choose 
to take an approach similar to that of semi-analytic models 
and take advantage of the sensitivity of the results to the 
wind parameters. By varying the parameters with halo mass 
or with the physical properties of the star- forming gas, one 
could match a wide range of SFHs. 

While we have not tried to tune the SFH, we have in- 
vestigated a toy model that uses the same amount of SN 
energy as the reference model, but in which the wind ve- 
locity scales with the local effective sound speed, Cs ,eff , as 
might be the case for thermally driven winds. If it is indeed 
hydrodynamic drag that stalls low-velocity winds, then this 
scaling could keep the winds efficient at all pressures. We 
implement this model, which we term WDENS, by making 
the wind parameters functions of the density of the gas from 
which the star particle formed: 



/ \i/6 




which implies oc Cs,cff since star-forming gas follows the 
effective EOS P = /OgC^^g cx Pg'^^ We set 7;* = 600 kms'^ 
and ry* = 2, so that the values of the wind parameters agree 
with those of the reference model for stars formed at the 
density threshold njj, while the wind velocity is greater (and 
the mass loading smaller) at higher pressures. 

Comparing the 25 Mpc runs shown in Fig. 1141 
we see that WDENS predicts a nearly identical SFH as 
WMLl V848 and REF down to z = 4, but that it generates 
much more efficient winds at later times. The 100 Mpc 
shows qualitatively the same behavior, with WDENS pre- 
dicting significantly lower SFRs below z = 2. Because the 
winds in WDENS remain effective for higher galaxy masses, 
the drop in the SFR below redshift 2 is steeper than for the 
reference model, although it is still less steep than observed. 

Comparing the two WDENS runs, we see that the 
agreement between the different box sizes is much worse 
than for the reference model. Clearly, making SN feedback 
a function of the local gas density increases the sensitivity 
to numerical resolution as we already concluded from the 
models that used a top-heavy IMF at high densities (see 
section l4.7.2|l . This is probably because the high density tail 
of the PDF can only be sampled if the galaxy contains a 
sufficient number of particles. One would therefore expect 
somewhat better convergence if the wind parameters were 
a function of the properties of the dark matter halo rather 
than the local gas properties. We will investigate such mod- 
els in Section [4.91 
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Figure 15. As Fig. \5\ but comparing the SFHs for models 
that use different implementations of SN feedback. Model WHY- 
DRODEC is identical to the reference model, except that the 
wind pa rticles are temporarily decoup led from the hydrodynam- 
ics as in lSpringel fc Hernguis^ ll2003al l. Model WTHERMAL, on 
the other hand, in jects the SN energy in therma l form following 
the prescription of lDalla Vecchia fc Schavel l l2009t) . Hydrodynam- 
ically decoupled wind particles can freely escape the ISM, but are 
unable to drag other particles along. They are therefore less effi- 
cient at high redshift, when low-mass galaxies dominate the SFR, 
but they are much more efficient at low redshift, when hydrody- 
namical drag within the high-pressure ISM of massive galaxies is 
important. Injecting the same amount of SN energy in thermal 
form increases the efficiency of the feedback in poorly resolved, 
low-mass galaxies, but the winds are somewhat less effective at 
low redshift, i.e. for higher mass galaxies. 

4-8.2 Hydrodynamically decoupled winds 

In recent years, a large fraction of the results from cosmo- 
logical, SFH simulations discussed in the literature wore ob- 
tained from simulat ions run with GADGET 2 (Springcl 200^ 
and employing the ISprineel fc HernauistI (|2003ah prescrip- 
tion for kinetic SN feedback. This prescription for galactic 
winds differs in two respects from ours. First, the wind parti- 
cles are selected stochastically from all the star-forming (i.e. 
dense) particles in the simulation and are therefore not local 
to the star particles as is the case for us. Second, the wind 
particles are subsequently decoupled from the hydrodynam- 
ics for 50 Myr (i.e. 31 kpc if traveling at 600 kms~^) or until 
their density has fallen below 10 per cent of the threshold 
for SF, which ensures that they esc ape the ISM. 

iDalla Vecchia fc Schavel l|2008l ) investigated the effects 
of this decoupling in detail and found them to be dramatic. 
While decoupled winds remove fuel for SF, they cannot blow 
bubbles in the disc, drive turbulence or create channels in 
gas with densities typical of the ISM. Decoupled winds are 
less efficient at suppressing SF in low mass galaxies, because 
they cannot drag gas along. They are, however, much more 
efficient for high mass galaxies, because they do not suffer 
the large energy losses due to drag in the high-pressure ISM. 
As the numerical resolution is decreased, the disc responsible 
for the drag disappears and the two prescriptions tend to 
converge. Hence, the decoupled winds are less sensitive to 
resolution, essentially because they only need to resolve the 
Jeans scales at the relatively low density for which the wind 
particles are recoupled. 
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Figure 16. As Fig. [5] but comparing the SFHs for mod- 
els that all employ tlie r mal S N feedback using the method of 
iDalla Vecchia &: Schayj | |2009|) . All models inject 40 per cent of 
the available SN energy per unit stellar mass, but they use differ- 
ent temperature jumps for the heated gas. Models WTHERMAL- 
LOGTxpy use temperature jumps log^Q AT = x.y (e.g. LOGT7p5 
implies AT = lO'^ ^ K) . Model WTHERMALLOGT7p5 assumed 
identical parameters as model WTHERMAL shown in Fig. 1151 
All simulations use a 12.5 Mpc box and 2 X 256^ particles. 
Enforcing a greater temperature jump, which implies a smaller 
heating probability, makes the feedback less efficient at high red- 
shift, but more efficient at low redshift. However, the feedback 
efficiency is less sensitive to the temperature jump than it is to 
the wind velocity for the case of kinetic SN feedback. 



Fig. [15] shows the SFH predicted by a model in which 
the wind particles were deco upled from the hydrodynamic s 
in the manner described in ISpringel fc Hernguistl (|2003al ) 
( WHYDRODEC) (note that the wind particles were, how- 
ever, still local to the newly formed star particles). Com- 
pared with the reference model, the decoupled winds are less 
efficient at high redshift, because the wind particles cannot 
drag other gas particles out of low-mass galaxies. We expect 
this difference to increase for higher resolution simulations, 
because they can resolve the discs of such galaxies better. 
Decoupled winds are, however, much more effective at lower 
redshifts when the galaxies dominating the SFR are better 
resolved and more massive, leading to large energy losses due 
to drag in the fSM for the reference model. Consequently, 
the peak in the SFR shifts from z ~ 2 for the reference run 
to 2 ~ 4 for the decoupled winds. 



4.8.3 Thermal SN feedback 

In the previous sections we found that the results of simu- 
lations employing kinetic SN feedback are sensitive to the 
parameters of the subgrid model, even if the total energy 
in the wind is kept constant. It is therefore of interest to 
consider alternatives to kinetic feedback. In particular, ther- 
mal feedback would seem a natural choice as it would allow 
the simulation itself to determine the properties of the wind 
such as the mass loading, velocity and geometry, based on 
the properties of the starburst. Unfortunately, current sim- 
ulations of cosmological volumes cannot resolve the energy- 
conserving phase in the evolution of SN remnants, causing 
any thermal energy input to be mostly radiated away before 



it can be converted into kinetic form. This is why alterna- 
tives such as kinetic feedback have been developed in the 
first place. 

As discussed in iDalla Vecchia fc Schavd l|2009l ). this 
overcooling problem is mainly caused by the fact that in 
naive implementations of thermal feedback the ratio of the 
heated mass to that of the star particle is too large, which 
means the temperature jump is too small and hence that 
the radiative cooling times are too short. The problem can 
thus be overcome if the SN energy produced by a star par- 
ticle is injected in a sufficiently small amount of mass, such 
that its cooling time becomes long compared to the time 
scale on which the loca l gas density can change in response 
to the energy injection. iDaUa Vecchia fc Schavj (2009) pro- 
pose a stocha stic method, whic h generalizes a prescription 
introduced bv lKav et all (|2003l ). that uses the temperature 
increase of the heated gas, AT, and the fraction of the 
SN energy that is injected, fth, as parameters. These pa- 
rameters then determine the probability for a gas element 
neighbouring a newly f o rmed star particle to be heated. 
iDalla Vecchia fc Schavd l|2009D demonstrated that as long 
as the temperature increase is greater than a value that de- 
pends weakly on the resolution and the gas density, the fac- 
tor by which the feedback suppresses SF is insensitive to the 
value of the temperature increase. 

We have performed a run, WTHERMAL, in the 
25 Mpc box using t h e ther mal feedback prescription of 
iDaUa Vecchia fc Schavd l|2009l ), setting AT = lO'''^ K and 
fth = 0.4. The latter value agrees with the one implied by 
the product -qv^ of the parameters of the kinetic feedback 
used in the reference run. Substitut ing these parameter val- 
ues in to the equations presented in lDalla Vecchia fc Schavd 
(200I), we estimate that for our resolution the thermal 
feedback will be efficient up to at least the density nu ~ 
2 X lO'^ cm"'', which exceeds our SF threshold by more than 
three orders of magnitude. Per star particle, the average 
number of gas particles that receive SN energy is about 0.54. 

Fig. [15] shows that the thermal feedback, when imple- 
mented in this manner, is indeed effective at suppressing 
the SFR. At high redshift it is more efficient than the ki- 
netic feedback used in the reference model, which suggests 
that it results in higher mass loading factors for poorly re- 
solved galaxies. For z < 5 the SFR is higher than in the 
reference run, but the difference is always less than 0.2 dex 
and remains constant below redshift 3.5. 

Fig. 1161 shows the SFHs for three different models that 
all employ thermal feedback using identical amounts of SN 
energy (/th ~ 0.4), but different temperature jumps. From 



top-to-bottom at redshift 2, the models use AT 



10' 



lO'^'^ and 10^'* K, respectively. Hence they differ by factors 
of two, which matches the factors of \/2 difference between 
the wind velocities used in Fig. 1141 These three thermal feed- 
back models use 2 x 256"' particles in a 12.5 Mpc box, 
which means the resolution is identical to that used in the 
L025N512 runs shown in Fig. [T4l 

Comparing the three models, we see that a higher tem- 
perature increase makes the feedback slightly less efficient 
at high redshift [z > 6), although the effect is marginal. 
Interestingly, the small difference at 2 > 6 appears to be 
be caused by the varying strength of the response to the 
reheating associated with reionization (which happens at 
Zr = 9 in our simulations). This agrees with the finding 
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of IPawlik fc Schavd (|2009l ') that SN feedback and photo- 
heating strengthen each other. Near the end of the simu- 
lations the feedback is more efficient for greater AT. The 
redshift at which AT — 10^ * K becomes more efficient than 
AT = 10^'^ K is lower than the redshift at which the latter 
model becomes more efficient than AT = lO''-^ K. All of this 
can be explained if the thermal feedback becomes inefficient 
for haloes more massive than some value which increases 
with AT. This situation parallels that of the kinetic feed- 
back with AT playing the role of «„• 

Interestingly, the SFH is much less sensitive to the value 
of AT than it is to «w for the case of kinetic SN feedback. 
The difference between the models with AT = 10^'^ and 
lO'^'* K is always less than 0.24 dex, whereas the differ- 
ence between the models with = 424 and 848 kms"'^ 
is 0.38 dex at z = 2 (see Fig. [14]). Moreover, while the SFHs 
in the different kinetic feedback models diverge rapidly, the 
differences between the thermal models is nearly constant 
below z — 3. 

Our findings that the thermal feedback is efficient and 
that it is less sensitive to the parameters of the model than is 
the case for kinetic feedback are very encouraging. We only 
used 40 per cent of the SN energy because we wanted to 
match the fraction used for the kinetic feedback in the ref- 
erence model. Higher values can, however, easily be justified 
for thermal feedback since we are now simulating radiative 
losses. 



4.9 "Momentum-driven" winds 

We have so far only considered galactic winds driven by 
SN feedback. It is, however, also possible that outflows 
are driven by the momentum that is deposited by photons 
from massive stars as they are absorbed by dust grains 
iMurrav et all |2005| ). In such a momentum-driven wind 
one would expect the mass loading to be inversely propor- 
tional to the wind v elocity and, according to the model of 
[Murray et all (|200d ). the terminal wind velocity w ould be 
simila r to the velocity dispersion of the host galaxy. iMartinI 
(|2005l ) argued that such scalings are in agreement with es- 
timates of the mass and velocit y in cold clouds as infe rred 
from Nal absorption lines. .Oppenheimer fc Pavel (|2006l ) im- 
plemented several versions of momentum-driven winds into 
cosmological simulations of the chemical enrichment of the 
IGM and found good agreement with observations of C IV 
absorption. 

Direct simulation of radiation pressure requires radia- 
tive t ransfer, which is too costly f or cosmological simula- 
tion s. lOppenheimer fc Pavel ll20q6f) therefore chose to use 
the ISpringel fc Hernguistl ( 2003al ') recipe for kinetic feed- 
back, including the hydrodynamic decoupling discussed in 
section l4.8.2l but to make the velocity kick, «w, and the mass 
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loading factor, r;, functions of the properties of the galaxy 
so as to mimic the scalings expected for momentum-driven 
winds. Speciflcally, they assumed 



Vv = ^ai-\/ a2f-L{Z) — 1 + as ) a, 



(7) 
(8) 



where ai — a4 are free parameters, a is the galaxy velocity 
dispersion, and fh{Z) is a function that accounts for the 
dependence of the stellar luminosity on metallicity and that 



.000 



^ 0.100 



0.010 



0.001 



; wpo 


REF - 

NOKICK ; 

WPOT 

WVCIRC 




a : 

\ \ - 
\ \ 

A \ 
\ \ ■ 



2 3 
Redshift 



5 6 7 8 9 



Figure 17. As Fig. [5] but comparing the SFHs for differ- 
ent im plementations of "momentum-driven" winds, following 
lOppen hcimcr & Dave ( 2006 , 2008) . While these models were mo- 
tivated by the claim that galact ic outflows may be driven by ra- 
diation pressure on dust grains iMurrav et al.ll2005f) , they do not 
actually include radiative transfer. Instead, the parameters of the 
kinetic feedback prescription are made functions of either the lo- 
cal potential {WPOT and WPOTNOKICK) or the mass of the 
dark matter halo ( WVGIRG). The wind parameters scale such 
that the wind velocity increases with halo mass, while the mass 
loading is inversely proportional to the wind velocity. For massive 
galaxies the winds use more energy than is available from SNe. 
See the text for additional details. Because the wind velocity in- 
creases with halo mass, the feedback is more efficient than for the 
reference model at late times. At high redshift it is also somewhat 
more eflicient thanks to the higher mass loading factors. 



varies from 1.7 at 10 ^ Zq to unity for solar abundances. 
Th e parameter a \ was set to 3 which is the value suggested 
bv iMurrav et al.l l|2005h . Parameter 02 was either set to 2 
or varied randomly between 1.05 and 2. Parameter 03 was 
introduced after noting that for radiatively driven winds the 
outflow velocity increases out to large distances, whereas in 
the simulations the gas is not given any more momentum 
after it is 'kicked' out of the ISM. They tried both a-j, = Q and 
as = 2. Finally, parameter 04 was set to 300 kms~^ in order 
to roughly match the observed cosmic SFR at high redshift. 
The velocity dispersion of the host galaxy was estimated by 
taking the local gravitational potential, "1>, and estimate a 
using the virial theorem {a = — 1$). 

Later papers by the same authors compared to other 
types of obse rvations, but used different parameter values 
and methods. lOppenheimer fc Pavel l|2008l ) realized that the 
gravitational potential is dominated by large-scale struc- 
ture rather than by the mass of individual haloes. For 
this reason, they moved away from using the gravita- 
tional potential to estimate wind properties and instead 
used friends-of-friends halo catalogues, generated on-the- 
fly throughout the simulation, and set a = \/2vc where 
Vc = \/GM / Rv ir is the halo circular veloci ty and Tivir is 
the virial radius. [Oppenheimer fc Pavel (I2OO8I ) also increased 
ai from 3 to 4.3, halved the value of 04 to 150 kms~^ 
and imposed an upper limit on the wind velocity corre- 
sponding to twice t he total SN energy. L a ter pa pers used 
the same values as lOppenheimer fc Pavel l|2008l ) although 
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lOppenheimer fc Pavel (|2009l ) no longer imposed a limit on 
the total wind energy. 

We first implemented the lOppenheimer fc Davi (|2006l ) 
method, which uses the local potential to estimate the ve- 
locity dispersion, but without decoupling the wind particles 

from the hydrodynamics. We used the parameter values ad- 

I 1 i — ^ 

vocated by lOppenheimer fc Davg 1120081 ) , although we ne- 
glected the metallicity-dependent square root term, which 
simplifies equation ^ to v„ = (ai + 03)0-. We ran two ver- 
sions. While both assumed ai = 3 and 04 — 150 kms^^, 
model WPOTNOKICK used as = whereas model WPOT 
used as = 2. Fig. llTl shows the SFHs predicted by the models 
for both box sizes. WPOT gives lower SFRs than WPOT- 
NOKICK which is not surprising since it uses higher wind 
velocities. Model WPOTNOKICK predicts higher SFRs 
than the reference model for the 25 Mpc box, but the 
order is reversed for the 100 Mpc box. This reflects the 
fact that there is more large-scale structure in the larger box 
and that the potential is dominated by the largest structures 
in the box. Clearly, this situation is not desirable. 

We therefore also ran model WVCIRC which estimates 
a using an on-the-fly ha l o frie nds-of-friends halo flnder as 
in lOppenheimer fc Pavel (|2008l f°l. The halo flnder was run 
at times spaced evenly in log a with Aa = 0.02a, where 
a is the expansion factor. Haloes were found based on the 
distribution of dark matter particles using a linking length 
of 0.2 and requiring a minimum of 25 dark matter particles 
per 

halcEl. 

Baryonic particles were attached to the nearest 
dark matter particle. 

The predicted SFHs are also shown in Fig. [T7| Except 
at very high redshift, the SFR is strongly reduced compared 
to the reference model. The difference increases with time, 
so that the SFR falls off more rapidly below z = 2 than for 
REF, as required by the observations. As was the case for 
WDENS (see Section g]!] and Fig. [HI), the increased effi- 
ciency of the feedback at late times arises because the wind 
velocity increases with the halo mass, whereas it is constant 
for the reference model. Contrary to WDENS, the feedback 
is also more efflcient than that of the reference model at 
high redshift. This is a consequence of the high mass load- 
ing used for low-mass haloes (recall that for WDENS the 
mass loading is never higher than for REF). 

The implementation of momentum-driven winds is 
rather crude. For example, the wind reaches its maximum 
velocity when it leaves the ISM rather than in the outer 
halo as expected for radiatively driven winds. Moreover, the 
parametrization leaves a lot of freedom, even more than for 
SN feedback, partly because the total amount of energy is 
no longer limitetf^. More to the point, Haas et al. (in prepa- 
ration) show that the models inject much more momentum 
than is actually available in the form of star light. It would 
therefore be dangerous to use comparisons between obser- 
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Note that lOppenheimer fc Pavel 1I2OO8I ) ran a friends-of-friends 
halo finder on the baryons using a linking length of 0.04 and as- 
suming a fixed baryon to dark matter ratio equal to the universal 
value. 

For star particles forming outside of any halo we assume the 
minimum halo mass corresponding to 25 dark matter particles. 

Depending on redshift, WVCIRC injects more energy in the 
winds per unit stellar mass formed than model REF for halo 
masses that exceed 10^^ — 10^^ M©. 
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Figure 18. As Fig. [Sj but comparing the SFHs for models with 
and without supermassive BHs. Feedback from AGN strongly 
suppresses the SFR in massive galaxies, which becomes more im- 
portant at late times. When AGN are included the SFR appears 
to be too low compared with observations, but that could be 
changed by decreasing the fraction of the SN energy that is in- 
jected, which was in this case chosen to roughly match the peak 
of the SFH in the absence of AGN feedback. 



vations and models such as these to discriminate between 
outflows driven by SNe and radiation pressure. Moreover, 
we note that it is not clear that observations of high-mass 
galaxies should agree with the "momentum-driven" models, 
given that AGN feedback is thought to be crucia l for such 
objects (e.g. lCroton et al.|[2006l : iBower et al.ll2006l i. but was 
not included. Conversely, if simulations such as the ones pre- 
sented here were to disagree with observations, it would not 
necessarily mean that winds are not radiatively driven. 

It is interesting that the scalings of the "momentum- 
driven" prescription, i.e. a wind velocity that increases and 
a mass loading that decreases with galaxy mass, agree quali- 
tatively with the results obtained fo r high-resolution simula- 
tions o f ener gy-driven SN feedback l|Dalla Vecchia fc Schavel 
|2008| . l2009| j. Even when hydrodynamical interactions 
are not temporarily ignored, the mass loading will be 
underestimated for galaxies that are poorly resolved 
dPalla Vecchia fc Schavel I2OO8I) . It may therefore be that 
the "momentum-driven" wind scaling partly compensates 
for some of the resolution effects that plague cosmological 
simulations, particularly at high redshift. 

4.10 AGN Feedback 

The centers of galaxies are thought to harbor supermassive 
BHs. Matter accreting onto these BHs emits large amounts 
of high energy radiation. Even if only a small fraction of this 
energy gets coupled to the ISM, it could have a dramatic ef- 
fect on the galaxies. Moreover, the magnetic flelds carried 
by the accreting matter could lead to the formation of jets 
which can displace and heat gas in and around galaxies. 



the low SFRs of high-mass galaxies (e.g. 


Pi Matteo et al. 


2005l:ICroton et al.ll2006l:lBower et al.ll2006l: 


Booth fc Schave 



I2009al ') and the suppression of cooling flows in clusters 
of galaxies (e.g. Churazov et al. 200 ll: Palla Vecchia et akl 



I2OO4I : iRuszkowski et al.ll2004l : ICattaneo fc Tevssiedl2007l ) 



© 2007 RAS, MNRAS OOO.nWl 



The physics driving the star formation history 23 



To investigate the effect of AGN feedback, we have run 
a series of simulations employing the subgrid prescription 
for the growth of BHs an d feedback from AGN described in 
iBooth fc Schavj (|2009al ) which is a substantia lly modified 
version of the model of ISpringel et al] (|2005al ). Below we 
will briefly summarise the main feature s of this model, but 
we refer the reader to iBooth fc Schavd (|2009al ) for further 
details and tests. 

Seed BHs of mass msecd are placed into every dark mat- 
ter halo whose mass exceeds mhaio.min- Our fiducial AGN 
model uses msocd = 9 x 10'' M0 , which corresponds to 
10"^ mg (6.4 X 10"^ mg) in the 100 h'^ Mpc (25 h'^ Mpc) 
box. For the minimum dark halo mass we use mhaio.min = 
4 X 10'° Mq, which corresponds to 10^ (6.4 x 10^) dark 
matter particles in the 100 Mpc (25 Mpc) box. 
Haloes are identified by running a friends-of-friends group 
finder on-the-fly as described in section 14.91 BHs can grow 
via Eddington-limited accretion of the surrounding gas and 
thr ough mergers with other BHs. 

iBooth fc Schavd (1200931 ) show that any simulation that 
resolves the Jeans scales will also resolve the Bondi-Hoyle 
accretion radius for BHs whose mass exceeds the simula- 
tion's mass resolution. At densities below the SF threshold 
rig (10~^ cm~^ in our model), the gas is kept warm by the 
lonizmg background (|Schavjf20o3 ) and we marginally re- 
solve the Jeans scales in our highest resolution runs (see 
Section [3.2p . For higher densities, however, a cold phase is 
expected to be present and naive application of the Bondi- 
Hoyle formula would lead us to strongly underestimate the 
accretion rate. We therefore assume that the accretion rate 
is given by the minimum of the Eddington rate and 
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47rG m^iip 
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where men is the mass of the BH, Cs and p are the sound 
speed and density of the local medium, v is the velocity of 
the BH relative to the ambient medium, and a is a dimen- 
sionless efficiency parameter given h^^^ 



( 1H 



1 if riH < 

otherwise. 



(10) 



Observe that for q = 1 (i.e. nn < nn or /9 = 0), equation ^ 
reduces to the Bondi-Hoyle-L yttleton rate l|Bondi fc Hovld 
1 19441 : iHovle fc Lvttletonlll939l ). Our fiducial AGN runs use 
/3 = 2, which results in efficient BH growth in haloes with 
stellar masses >1O'°-^M0 in the LlOO runs. 

The amount of accreted mass is related to the rate 
of growth of the BH by riiBH = nT-accr(l — tr), where er 
is the radiative efficiency of a BH, which we always as- 
sume to be 10%, the mean value for the radiatively efficient 
IShakura fc SunvaevI (|l973l ) accretion onto a Schwarzschild 
BH. 

We assume that a fraction et of the radiated energy 
couples to the ISM. The amount of energy returned by a 
BH to its surrounding medium is thus given by 



" Note that ISpringel et"al] ll2005a') used a fixed value a = 100. 
Consequently, massive BHs need to suppress the ambient gas den- 
sity to values far below i n order to reduce the ac cretion rate 
to sub- Eddington values fsee lBooth &: Schavd[2009a^ . 



where c is the speed of light. We set et = 0.15 in order to 
match the observed cosmic mass density in BHs as well as 
the relation between BH and galaxy mass, both at redshift 
zero. BH particles store feedback energy until it suffices to 
heat Tihcat of their neighbours by ATmin. The two param- 
eters ATmin and Wheat are chosen such that AGN heated 
gas obtains a long cooling time and so that the time taken 
to perform a feedback event is shorter than the Salpeter 
time for Eddington-limited accretion. The parameter choices 
ATinin = 10* K and rihcat = 1 are found to provide a good 
balance between these two constraints. 

Fig. [18] shows that the addition of AGN feedback 
strongly suppresses the SFR. The difference with the ref- 
erence model increases with time, which implies that AGN 
feedback is more important for higher mass galaxies. The 
drop in the SFR below z = 2 is much closer to the observed 
slope when AGN feedback is included, but the overall am- 
plitude of the SFR is probably too low. It is, however, im- 
portant to note that the observed SFRs are subject to large 
systematic uncertainties. Indeed, McCarthy et al. (in prepa- 
ration) find that the predicted stellar masses are in fact in 
good agreement with observations of groups of galaxies. 

Even if the simulation with AGN feedback really did 
form too few stars, it would not be a concern here. As dis- 
cussed in section 13.1.51 the fraction of the SN energy that 
is injected was chosen to roughly match the peak in the ob- 
served SFH. It is therefore unavoidable that including an 
extra form of efficient feedback, without making any other 
changes, reduces the SFR to values that are lower than ob- 
served. It would have been possible to adjust the parameters 
of the prescription for SN feedback to obtain a better match 
to t he observed SFH , but th is is not the objective here. 

iBooth fc Schavd l|2009al ) have carried out an extensive 
parameter study of the AGN model at the resolution of 
our 100 Mpc box, comparing the cosmic SFH as well 
as other observables. They found that the results are sensi- 
tive to the accretion model, i.e. the value of /3, which sets the 
halo mass above which the BHs can grow onto the scaling 
relations. Remarkably, they found that the SFH is nearly 
completely ind ependent of the feedbac k efficiency, but, in 
agreement with lDi Matteo et al.l l|2005h . the BH masses are 
inversely proportional to et- They explained this in terms of 
self-regulation: the BHs grow until they have injected suf- 
ficient energy to balance the infall driven by gas accretion 
onto haloes and radiative cooling (see also Di Matteo et al.l 
I2OO5I : [Robertson et al.|[2006l : iBooth fc Schavd[2009bh . If haff 
as much energy is injected per unit accreted mass, then the 
BHs need to grow twice as massive in order to inject the 
same amount of energy. Because the factor by which SF is 
suppressed depends on the amount of energy that is injected 
by the BHs, the SFH is insensitive to variations in the as- 
sumed efficiency of AGN feedback. 



5 SUMMARY AND DISCUSSION 

The cosmic star formation history (SFH) is perhaps the most 
fundamental observable in astrophysical cosmology. It is dif- 
ficult to model, because of the large range of galaxy masses 
that contribute and because of the many feedback processes 
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that may be important. Cosmological hydrodynamical simu- 
lations need to resort to subgrid prescriptions for the physics 
that remains unresolved, which makes them resemble semi- 
analytic models in some respects. It is, however, much more 
difficult to explore parameter space using fully numerical 
simulations because of the high computational cost. 

Here, we have introduced the Overwhelmingly Large 
Simulations (OWLS) project, which consists of more than 
50 large, cosmological, hydrodynamical simulations. The 
simulations were all run with a modified version of th e 
SPH code GADGETS (last described in iSpringell l2005ll. 
using new modules for radiative cooling ( Wiersma et al.l 
l2009al). SF JS chave fc Dal la Vecch"ial |2008|). chemodynam- 
ics |wiersma et al. 2009b'l .\inetic jPalla Vecchia fc Schavd 
I2OO8I ) or thermal (|Dalla Vecchia fc Schavd |2009| ) SN feed- 



back, and accretion onto and feedback from super massive 
BHs (jSpringel et al ] l2005al : iBooth fc Schav3 l2009al ) . With 
2 X 512'^ particles, the OWLS runs are among the largest 
dissipative simulations ever performed. The simulations are 
repeated many times, each time changing a single aspect of 
the input physics or a single numerical parameter with re- 
spect to the reference model described in section[3l We stress 
that this model merely functions as a reference point for our 
systematic exploration of parameter space and should there- 
fore not be regarded as our "best" model. 

Generically, we find that SF is limited by the build- 
up of dark matter haloes at high redshift, reaches a broad 
maximum at intermediate redshift, then decreases as it is 
quenched by lower cooling rates in hotter and lower density 
gas, gas exhaustion, and self-regulated feedback from stars 
and black holes, in broad agreement with previou s work (e.g. 
IWhite fc FrenMll99ll : iHernguist fc Springelll20o3 ) . While we 
have limited ourselves to comparisons of the cosmic SFHs, 
future papers will investigate the behavior of the models in 
other contexts. 

The models are run down to redshift 2: = 2 in boxes 
of 25 comoving Mpc on a side and/or down to 2: = 
in boxes of 100 Mpc. The mass (spatial) resolution is a 
factor 64 (4) better in the smaller box. Detailed convergence 
tests showed that a 25 Mpc box is sufficiently large to 
model the cosmic SFH down to 2 = 0. For the 25 /i"^ Mpc 
(100 Mpc) simulation of the reference model, the SFH 
is very close to converged for 2 < 7 (2 < 3). Interestingly, 
we found that while the star formation rate (SFR) typically 
increases if the resolution is improved, the situation reverses 
at low redshift before convergence has been attained. This 
calls into question the strategy to combine different box sizes 
to obtain a converged SFH all the way from high redshift to 

2 = 0. 

The cosmic SFR density can be decomposed into a halo 
mass function and the (distribution of the) SFR as a func- 
tion of halo mass. The baryonic physics mainly affects the 
latter function, which we shall explore in more detail in Haas 
et al. (in preparation) , while the mass function is determined 
mostly by the assumed cosmology. We showed that the SFH 
is sensitive to even relatively small changes in the cosmolog- 
ical parameters, such as the difference between the values 
inferred from the WMAP 3-year data used here and the ear- 
lier values assumed in the Millennium simulation. Clearly, 
the parameters of semi-analytic models build onto the Mil- 
lennium simulation will have to be modified if they are used 
on a simulation assuming the current concordance cosmol- 



ogy. Comparisons to observations of rare objects and the 
high-redshift Universe will be particularly strongly affected 
due to the relatively large difference in the value of og,- 

Our systematic tests of the subgrid physics revealed 
that SF in intermediate mass galaxies is highly self-regulated 
by feedback from massive stars. This explains our remark- 
able finding that the predicted SFH is nearly completely in- 
dependent of the treatment of the unresolved ISM, including 
the assumed SF law. If the SF efficiency is increased, then 
the galaxies simply reduce their gas fractions so as to keep 
the SFR, and thus the rate at which stars inject energy into 
the ISM, constant. Similarly, if the efficiency of SN feedback 
is changed by injecting a different fraction of the SN energy, 
then, to first order, the galaxies simply adjust their SFRs so 
as to keep the rate of energy injection constant. The critical 
rate of energy injection that results from self-regulation is 
presumably the rate required to balance gas infall resulting 
from accretion onto haloes and radiative cooling and will 
therefore depend on the halo mass and redshift. 

Note, however, that our findings apply to SFRs that 
have effectively been averaged over entire galaxies and over 
very long time scales. Self-regulation may also occur on 
smaller length and time scales. Indeed, we implicitly as- 
sume this to be the case through our use of empirical 
SF laws that have been averaged over spatial scales that 
are large compared to individual star-forming regions. Pro- 
cesses other than SN feedback may be important for the 
self-regulation that occurs on these smaller scales and per- 
haps even on large scales if SN feedback is ineffi cient (e.g. 
iRobertson fc Kravtsovll2008l : iGnedin et al.ll2009l ). 

For massive galaxies feedback from massive stars be- 
comes inefficient and it is the self-regulated growth of su- 
permassive BHs that ensures that a fixed amount of en- 
ergy is inj ected into the ISM (e.g. iDi Matteo et al.1 l2005l: 
iBooth fc Schavc 2009a b). As we showed in lBooth fc Schayg 
(2 009a !), the SFH is therefore independent of the assumed 
efficiency of AGN feedback. If the BHs inject twice as much 
energy per unit accreted mass, then they grow only half as 
massive so that they inject the same amount of energy and 
because they inject the same amount of energy, the SFR is 
suppressed by the same factor. 

At very high redshift the SFH is not yet controlled by 
self-regulation and will thus be limited by the gas consump- 
tion time scale implied by the assumed SF law. In our simu- 
lations the dependence on the SF law at high redshift may, 
however, mostly result from our limited resolution. Stars can 
only form in haloes that exceed the resolution limit and the 
formation of the first generation of stars within a halo is not 
hampered by winds. Because the minimum halo mass is set 
by the resolution, so is the SFH until it is dominated by 
haloes that are significantly more massive. 

Radiative feedback from non-local SF is, however, very 
important at high redshift. In particular, reheating associ- 
ated with hydrogen reionization quenches SF in haloes with 
virial temperatures < 10'* K and thus strongly affects the 
SFH when it would otherwise be dominated by such low- 
mass haloes. Our simulations underestimate this effect due 
to their limited resolution and because we assume that a 
photo-dissociating background is present at all redshifts. 
Contrary to hydrogen reionization, the milder reheating as- 
sociated with helium reionization does not have a significant 
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impact on the SFH because it is limited to the low-density 
IGM. 

While cosmology and self-regulation via outflows are the 
principal ingredients controlling the SFH, there are other 
processes that are important. Metal-line cooling becomes 
very important at late times, as the metal content of the 
gas builds up and less gas falls in cold. The inclusion of 
cooling by heavy elements has two important effects. First, 
it allows more virialized gas to cool, an effect that has been 
widely discussed in the literature, ffowever, we find that the 
dominant effect of metal-line cooling may be that it reduces 
the efhciency of galactic winds, as it increases thermal losses 
in the gas that has been shock-heated by the winds. The 
inclusion of metal-line cooling makes it much more difRcult 
to reproduce the steep decline in the observed SFH below 

2 = 2. 

Another process that becomes important at late times, 
is mass loss by intermediate mass stars. On time scales of 
hundreds of millions to billions of years, winds from AGB 
stars boost the SFR by providing fresh fuel for SF and by 
releasing metals that were locked up in stars into the ISM. 
As was the case for metal-line cooling, including mass loss by 
AGB stars make it more difficult to match the sharp drop in 
the cosmic SFR that is observed at low redshift. The shape 
of the time delay function for SNe of type la, on the other 
hand, turned out to be unimportant, presumably because 
its impact is limited to a shift in the timing of the release of 
a fraction of the iron. 

Without AGN feedback, it is challenging to match the 
steep decline in the cosmic SFR below 2 = 2. It is, however, 
not yet totally clear that feedback from massive stars cannot 
solve the problem. If injected in kinetic form, SN feedback 
becomes inefficient once the wind velocity falls below some 
critical value that increa ses with galaxy mass and thus with 
the pressure in the ISM l|DaUa Vecchia fc Scha"v3l2008h . By 
increasing the input wind velocity with galaxy mass, or with 
some correlated local property (e.g. gas density, pressure, or 
velocity dispersion) , while decreasing the input mass loading 
so as to keep the wind energy constant, we can keep SN 
feedback efficient in higher mass galaxies. 

Using a sufficiently high constant wind velocity also 
results in efficient suppression of SF in relatively massive 
galaxies, but not for poorly resolved low-mass galaxies. If a 
galaxy is poorly resolved, then there is no disc from which 
wind particles can drag gas along, li miting the effective mass 
loading factor to the input value l|Dalla Vecchia Sz Schavel 
I2OO8I ) . Decreasing the input wind velocity (and thus increas- 
ing the input mass loading) with decreasing galaxy mass 
counteracts the resulting underestimate of the mass load- 
ing factor. If we allow ourselves to vary the parameters of 
the kinetic feedback with local properties, then we can even 
"design" SFHs, but this freedom comes at the expense of in- 
troducing additional free parameters and an increased sen- 
sitivity to the numerical resolution. 

Recently, the possibility that galactic outffows are 
driven by radiation pressu re on dust grains has gener- 
ated considerable interest (|Murrav et al.l l2005h . For such 
winds the mass loading is expected to be inversely pro- 
portional to the wind velocity, which, in the outer halo, 
is ex pected to be similar to the gala xy velocity disper- 
sion. lOppenheimer fc Pavel l|2006l , I2OO8I ) implemented such 
"momentum-driven" winds in a simplified fashion by kick- 



ing particles out of the ISM with several times the velocity 
that the wind is expected to reach in the outer halo, and 
by tuning the normalization of the initial wind mass load- 
ing to match the observed SFR. We ran several simulations 
that employed such methods and found that the feedback is 
more efficient than our standard SN feedback, particularly 
at low redshift when relatively massive galaxies dominate 
the SFR. This is partly because the winds are allowed to 
carry more energy than is available from SNe, but it may be 
mostly due to the fact that the "momentum-driven" scalings 
happen to overcome some of the numerical effects discussed 
above. Given the simplified nature and the limitations of the 
kinetic feedback models, it would be dangerous to use them 
to discriminate between winds driven by SNe and radiation 
pressure. Finally, we note that Haas et al. (in preparation) 
demonstrate that the amount of momentum that is in jected 
in the models of lOppenheimer fc Pavel lj2006l . I2OO8I ). and 
thus also in our versions of these models, exceeds the amount 
of momentum that is actually available in the form of star 
light by up to an order of magnitude. 

Besides the values of the wind parameters, the imple- 
mentation of kinetic SN feedback is also important. For ex- 
ample, most studies in the literature employing GADGET 
temporarily decouple wind particles from the hydrodynam- 
ics, allowing them to freely escape the ISM without blowing 
bubbles or generating turbulence. We found that this de- 
coupling (which we included in one of our runs) , reduces the 
efficiency of SN feedback at high redshift, because of the in- 
ability of the winds to drag gas along, but strongly increases 
it at low redshift when hydrodynamic drag would otherwise 
quench the winds in massive galaxies. 

Outflows driven by feedback from SF, be it SN or radia- 
tion pressure, are very important. It is therefore unfortunate 
that the predictions for the SFH of simulations that imple- 
ment outflows in the form of kinetic feedback, i.e. by kicking 
particles, are not robust at the currently attainable resolu- 
tion. The predicted SFH is sensitive to the values of poorly 
constrained wind parameters, even for a constant wind en- 
ergy, and to the details of the implementation. Note that the 
same may well be true for other types of subgrid prescrip- 
tions for feedback from SF. Clearly, it is also crucial to vary 
the parameters of such models. On the other hand, we will 
show in future papers that many observables are much less 
sensitive to the implementation of feedback from SF than is 
the case for the cosmic SFH. 

We investigated implementing SN feedback in thermal 
form, using the method of P alia Vecchia fc Scha^ l|2009D to 
overcome the overcooling problem that is commonly encoun- 
tered when thermal SN feedback is used. Encouragingly, we 
found thermal feedback to be efficient and the predictions 
to be less sensitive to the values of the free parameters than 
is the case for kinetic feedback. Relative to kinetic feedback, 
thermal feedback is particularly efficient at high redshift. Be- 
cause heated particles push all their neighbours, the winds 
remain highly mass loaded even in poorly resolved galaxies. 

Another large uncertainty is the stellar IMF. The con- 
sequences of a change in the assumed IMF are difficult to 
predict because the IMF affects many things. The IMF de- 
termines the effective nucleosynthetic yields, the fraction of 
the stellar mass that is recycled by intermediate mass stars, 
and the amount of SN energy (and radiation) produced by 
massive stars. Moreover, a change in the IMF also implies 
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a change in the observed SF law and history, as SFRs are 
inferred from observations of the hght emitted by massive 
stars. These latter two effects are often ignored. Interest- 
ingly, we found that an IMF that becomes top-heavy at 
high gas pressures would improve agreement with the ob- 
servations because it preferentially suppresses SF in massive 
galaxies. 

How do we proceed from here? Based on our findings, a 
better understanding of the mechanisms responsible for the 
generation of galactic winds, as well as more robust numeri- 
cal implementations, are crucial. Paradoxically, a better un- 
derstanding of SF may not directly improve predictions for 
the cosmic SFH if, as is the case for the models studied here, 
galaxies regulate their SFRs by adjusting their gas fractions 
via large-scale outflows. At very high redshift photo-heating 
is key and simulations including radiative transfer would 
constitute a clear improvement. At lower redshifts it is very 
important to include mass loss from AGB stars and metal- 
line cooling. Better treatments of metal mixing and the in- 
clusion of non-equilibrium cooling would certainly be helpful 
here. The steep drop in the observed SFR below 2 = 2 is dif- 
ficult to reproduce without AGN feedback. Fortunately, the 
self-regulatory nature of the growth of supermassive BHs 
makes predictions for the SFH insensitive to the assumed 
efficiency of AGN feedback. Finally, as always, higher reso- 
lution would be very helpful as it would, for example, enable 
us to probe further down the halo mass function and to in- 
clude more realistic treatments of galactic winds. 
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